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FOREWORD 


NASTRAN® (NASA STRUCTURAL ANALYSIS) is a large, comprehensive, nonproprie- 
tary, general purpose finite element computer code for structural analysis 
which was developed under NASA sponsorship and became available to the public 
in late 1970. It can be obtained through COSMIC (Computer Software Management 
and Information Center), Athens, Georgia, and is widely used by NASA, other 
government agencies, and industry. 

NASA currently provides continuing maintenance of NASTRAN® through COSMIC. 
Because of the widespread interest in NASTRAN®, and finite element methods in 
general, the Ninth NASTRAN® Users' Colloquium was organized and held at the 
Kennedy Space Center, October 22-23, 1980. (Papers from previous colloquia 
held in 1971, 1972, 1973, 1975, 1976, 1977, 1978, and 1979 are published in 
NASA Technical Memorandums X-2378, X-2637, X-2893, X-3278, X-3428, and NASA 
Conference Publications 2018, 2062, and 2131.) The Ninth Colloquium provides 
some comprehensive general papers on the application of finite element methods 
in engineering, comparisons with other approaches, unique applications, pre- 
and post-processing or auxiliary programs, and new methods of analysis with 
NASTRAN. 

Individuals actively engaged in the use of finite elements or NASTRAN® 
were invited to prepare papers for presentation at the Colloquium. These 
papers are included in this volume. No editorial review was provided by NASA 
or COSMIC, however, detailed instructions were provided each author to achieve 
reasonably consistent paper format and content. The opinions and data pre- 
sented are the sole responsibility of the authors and their respective organi- 
zations. 

Cochairmen: 

Robert L. Brugh 

COSMIC 

University of Georgia 

Athens, GA 30602 

and 

Henry Harris 

John F. Kennedy Space Center 

Kennedy Space Center, FL 32899 
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NEW CAPABILITIES AND MODIFICATIONS FOR NASTRAN 


LEVEL 17.5 AT DTNSRDC 
Myles M. Hurwitz 

David W. Taylor Naval Ship Research and Development Center 


SUMMARY 


Since. 1970 DTNSRDC has been modifying NASTRAN to suit various Navy require- 
ments. These modifications have included major new capabilities as well as 
user conveniences and error corrections. This paper describes the new features 
added to NASTRAN Level 17.5 at DTNSRDC. The subject areas of the additions 
include magnetostatics, piezoelectricity, fluid-structure interactions, iso- 
parametric finite elements, and shock design for shipboard equipment. 


INTRODUCTION 


The David W. Taylor Naval Ship Research and Development Center (DTNSRDC) 
has been involved with NASTRAN since 1968. In the ensuing 3-4 years, prior to 
the first public release of the program in 1972, engineers at DTNSRDC gained 
valuable experience with NASTRAN, often interacting with the program developers 
on various theoretical, programming, and user aspects. The result of that 
early effort was a detailed NASTRAN evaluation report, which included a brief 
description of our first modification to NASTRAN- — the addition of a heat trans- 
fer finite element to the NASTRAN element library (ref . 1) . 

In subsequent years, the DTNSRDC modifications to NASTRAN were many and 
varied, ranging from error correction and user conveniences to new finite 
elements and new functional modules and rigid formats. 

Since Level 17.5 was released in the Spring of 1979, our NASTRAN modifica- 
tion effort has remained vigorous. The subject areas of new capabilities and 
uses include magnetostatics, piezoelectricity, fluid-structure interactions, 
isoparametric finite elements, and shock design of shipboard equipment. This 
paper describes these subject areas as we have implemented them into NASTRAN, 
sample applications of some of these areas, and the correction of an important 
program error. All of this work will be transferred to the DTNSRDC version of 
Level 17.6 after the standard version is released. 


MAGNETOSTATICS 


The prediction of static magnetic fields around ships 


and submarines is of 


1 


concern to the Navy because of the need for these craft to remain undetected. 

A numerical procedure which can predict these fields can also be used to 
evaluate systems which might reduce the fields, e.g. degaussing coils. Such a 
procedure, making use of a scalar potential rather than the less efficient 
vector potential, was described in reference 2. Reference 3 describes a 
capability for computing the magnetostatic fields about axisymmetric structures 
that was added to NASTRAN. However, that work was limited to the TRAPRG and 
TRIARG finite elements and to axisymmetric current coils. In Level 17.5, the 
analysis has been extended to general built-up and continuum structures with 
general current coil configurations. The finite elements allowed are those 
available for NASTRAN heat transfer analysis (ref. 4), for reasons which may be 
seen from the brief description of the applicable theory which follows. 

The applicable Maxwell equations governing the magnetostatic case are 

V xH = J (1) 


7 • B = 0 


( 2 ) 


where 

H = magnetic field strength or intensity 
B = magnetic induction or flux density 
J = current density 

The constitutive relation 


B = pH (3) 

where p is the magnetic permeability is also required. If H is separated into 
two parts 

H = H + H (4) 

c m 

where H , the field in a homogeneous region due to current density J (as might 
occur in a current coil), may be computed using the Biot-Savart law, (ref. 5), 
then H^j becomes the only unknown. By equations (1) and (4), 

V xH = 0 (5) 

m 

so that 

H = V<f) (6) 

m 

where <j> is a scalar potential. By equations (2), (3), (4), and (6), 

V • pVij) + V • pH^ = 0 (7) 

which can be put into the standard form 

K<j> = F (8) 
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where 


k. . = / (VN.) T p(VN.)dV (9) 

1 J v 1 J 

f. = -/ (vn.) t m dV (10) 

i v i c 

N. being the displacement function for a finite element at the i^ grid point. 
Equation (9) is of the same form as that required to compute the conductivity 
matrix in heat transfer, with p representing magnetic permeability rather than 
thermal conductivity. Equation (10) , which is dependent on the finite element 
type, is not in a standard heat transfer form and was added to NASTRAN along 
with the new bulk data cards needed to specify H c . Current coils may be 
defined, from which NASTRAN computes H c using the Biot-Savart law, or H c may be 
specified as coming from an ambient field, or a combination of both sources of 
H c may be given. 

Equation (6) gives the unknown Hjjj, which, in standard NASTRAN terminology, 
is the thermal gradient, and equations (4) and (3) yield the final result. 

One unanticipated addition to NASTRAN was required when it was discovered 
that the program did not compute thermal gradients for the isoparametric 
solids IHEX1 , IHEX2 , and IHEX3, as needed by equation (6). An example of this 
capability is shown in figure 1. The finite element model depicts a solid 
sphere (shaded part) which has been placed into a uniform, ambient, axial 
magnetic field. TRIARG elements were used and only the upper half was modeled 
due to symmetry. The NASTRAN results are compared with theoretical results in 
Table 1. 


PIEZOELECTRICITY 


The analysis of sonar transducers requires accounting for the effects of 
their piezoelectric materials. The theory for these materials couples the 
structural displacements to electric potentials (ref. 6). This theory was 
incorporated into NASTRAN only for the TRIAAX and TRAPAX finite elements (ref. 
7). These elements, trapezoidal and triangular in cross-section respectively, 
are solid, axisymmetric rings whose degrees-of-f reedom are expanded into 
Fourier series, thus allowing nonaxi symmetric loads. 

The piezoelectric constitutive relations may be written as 


'{a}\ 

_ 

•U E ] 

‘*>1! 

f{e> 

{D}j 

,U] T 

-u s J 1 

({E} 

where (a) = stress components = 

la , 
L rr 

a zz’ °ee’ 

°rz ’ °r0 * a zeJ T 
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{D} = components of electric flux density = [d , 


zz 


D 


60 


I 


{e} = mechanical strain components 

{E} = electric field components 

E 

[c ] = elastic stiffness tensor evaluated at constant electric field 

[e] = piezoelectric tensor 

g 

[e ] = dielectric tensor evaluated at constant mechanical strain 


The displacement vector of a point within an element is taken to be 



( 12 ) 


where u, v, and w are the ring displacements in the radial, tangential, and 
axial directions, respectively, and <J> is the electric potential. The latter 
degree-of-freedom is taken to be the fourth degree-of-f reedom at each ring. 
Each of these quantities is expanded into a Fourier series with respect to the 
azimuth position 0. The Fourier series for the electric potential has the 
same form as the Fourier series for radial displacement u, as given in the 
NASTRAN Theoretical Manual (ref . 4) . 

The "stiffness" matrix for the N c harmonic is 


[K <N) ] - „ JJ [b w ] 

r z 


(H),T f[c E ! [e] " 
_[e] T ~[e S ]_ 


[B (N) ] rdrdz 


. (N) 


where [Is ] is the matrix of "strain-displacement" coefficients for the N 
harmonic. 


(13) 

th 


Equations (12) and (13) indicate that the matrix equation to be solved for 
static analysis may be partitioned as follows: 
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{ V 

= vector of 

structural forces 



and 

{ V 

= vector of 

electrical charges 
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In addition to the new data cards describing the piezoelectric materials, 
many modifications and corrections to NASTRAN were made, including the computa- 
tion of complex stresses and forces for the TRAPAX and TRIAAX elements. 

An example of a piezoelectric problem is shown in figure 2. This is an 
axially polarized PZT-4 piezoelectric disk, whose natural frequencies are to be 
determined. Table 2 compares the NASTRAN results with experimental and MARTSAM 
results (ref. 8). MARTSAM uses finite elements similar to NASTRAN' s TRIAAX and 
TRAPAX elements, but with quadratic displacement functions rather than the 
linear displacement functions in NASTRAN. The MARTSAM results were obtained 
with a much coarser mesh. 


FLUID- STRUCTURE INTERACTION 


Investigation of . the coupling of fluid and structural effects has been an 
important part of the DTNSRDC program during the past few years. Applications 
include vibrations of submerged structures (refs. 9 and 10), shock response of 
submerged structures (refs. 11 and 12), and the response of fluid-filled pipes. 

Although these new applications did not require additions to NASTRAN, they 
did involve inventive use of DMAP and unusual use of existing data cards. This 
relatively new subject area shows the power of NASTRAN and its DMAP capability 
to adapt to new uses without requiring modification of the source code. 


ISOPARAMETRIC FINITE ELEMENTS 


A number of additions and modifications have been made to NASTRAN in the 
area of isoparametric finite elements. 

1. A two-dimensional membrane element IS2D8, with quadratic displacement 
functions, was added to the finite element library. This element has complete 
NASTRAN capability with the exception of piecewise-linear analysis. The 
element has been used in a number of applications where the CQDMEM1 element 
would have required a much finer mesh. 

2. The standard version of NASTRAN computes grid point stresses of the 
isoparametric solids IHEX1, IHEX2, and IHEX3 directly at the grid points. 
However, it has been shown that the stresses computed at the grid points are 
inferior to stresses extrapolated to the grid points from stresses calculated 
at the Gauss integration points (ref. 13). This extrapolation method has been 
added to the program for the IHEX1, IHEX2, IHEX3, and IS2D8 elements. 

3. The isoparametric solid elements are limited to isotropic materials in 
the standard version of NASTRAN. We have added a capability for rectangular 
anisotropy for those elements. 
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4. As mentioned in the Magnetostatics section, a thermal gradient 
computation has been added for the isoparametric solids. 

5. Although Level 17.5 allows for the choice of single precision or 
double precision arithmetic for some computations, including element matrix 
generation, it did not allow this choice for the isoparametric solids; only 
double precision was allowed. Since DTNSRDC uses CDC computers with 60-bit 
single precision words, a single-double choice for these elements was added. 
Generation time for the single precision stiffness matrix for one IHEX2 element 
with three Gauss integration points was reduced on the CDC 6400 computer from 
12 CPU seconds to 4 . 


SHOCK DESIGN OF SHIPBOARD EQUIPMENT 


The Dynamic Design- Analysis Method (DDAM) was developed for the shock 
design of shipboard equipment (ref. 14). This method is similar in many 
respects to the techniques used in earthquake analysis. In fact, an earthquake 
analysis using NASTRAN has been performed (ref. 15). However, DDAM, rather 
than some variation of it, is required by naval shipbuilding specifications for 
shipboard equipment. Therefore, we are presently developing a DMAP procedure 
and a functional module to perform DDAM analyses. 

Briefly, the steps in the DDAM method are as follows: 

1. Compute the normal modes and natural frequencies. 

2. For each mode, the i^, for example, compute the participation factor 

P. = — {<j>. } T [M]{D) (15) 

i m. i 
i 

= participation factor for the i*-* 1 mode 

X 

= modal mass term for the i^ mode = [M] { } 

= mass matrix 
= i^h mode shape 

= direction cosine vector defining desired direction (DDAM analyzes 
one direction at a time) 

3. Calculate the effective mass and effective weight in each mode. 


M e££ 

X 

= p.4. } T [M]{D} = m.P? 

X X XX 

(16) 

w e£f 

X 

M ef f 
= g M ± 

(17) 


where 


where P . 

x 

m. 

l 

[M] 

4 ± > 

(D> 
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= effective mass in i^ mode 

Gf f , « 

W = effective weight in i 4 - 11 mode 

g = acceleration of gravity 


4. Using the effective weights just computed, locate the design spectrum 
value for each mode in the desired direction. 

5. Compute the effective static force for each mode. 


{F.} = P.V.u. [M] {<* . } (18) 

x 1 i r 1 

th 

where uk is the i natural frequency. 

6. Perform a static analysis with each load to compute stresses. (There will 
be one static analysis for each desired mode in each desired direction.) 


7. Compute the so-called NRL sum (ref. 16) of the stresses at each desired 
point (element centroids) as follows: 


S. 

1 


S. 



(S.. ) 

Jb 


2 


where 


S. 


= the maximum stress at the jth point (taken over the modes under 
consideration) 


Two NASTRAN runs will be required for a complete DDAM analysis; the first 
will perform steps 1-3, and the second, steps 5-7. The D and V terms will be 
input through DMI cards, although some default values will be available for V. 

A post-processor, possibly included in NASTRAN as a new functional module, will 
perform the NRL sums in step 7. 


ERROR CORRECTIONS 

Numerous error corrections have been made to Level 17.5 by DTNSRDC and 
.reported to COSMIC, but perhaps the most important involved the stiffness 
matrix computation for the six elements QDMEM1, QDMEM2, SHEAR, TWIST, TRAP AX, 
and TRIAAX. The method of matrix computation for these elements was changed 
from SMA-type in Level 17.0 to EMG-type in Level 17.5. The change introduced 
an error which occurred only for certain combinations of grid point numberings 
for these elements . 

All the error corrections reported to COSMIC are expected to appear in the 
forthcoming Level 17.6. 
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I Table 1. Ferromagnetic Sphere Results 


- DATA ANALYTIC NASTRAN DEVIATION 

POINT SOLUTION SOLUTION 


1 

0.396 TESLA 56.3° 

0.396 TESLA 59.0° 

0.0% 

2.7 

2 

0.843 

52.7 

0.840 

53.5 

0.3 

0.8 

3 

1.488 

0.0 

1.537 

0.0 

3.3 

0.0 

4 

0.523 

86.9 

0.527 

87.1 

0.7 

0.2 

5 

0.941 

76.5 

0.921 

76.3 

2.1 

0.2 

6 

0.417 

79.2 

0.409 

80.6 

1.9 

1.4 

7 

0.571 

74.6 

0.579 

74.5 

1.5 

0.1 

8 

0.705 

83.7 

0.697 

83.6 

1.1 

0.1 

9 

0.479 

86.7 

0.483 

86.8 

.7 

0.1 

10 

0.526 

85.3 

0.532 

85.4 

1.1 

0.1 

11 

0.566 

88.3 

0.562 

88.5 

0.7 

0.2 

12 

0.492 

89.3 

0.499 

89.3 

1.5 

0.1 

13 , 

0.499 

88.6 

0.505 

88.6 

1.3 

0.0 

14 

0.516 

89.1 

0.523 

89.0 

1.6 

0.1 


Table 2. Natural Frequencies of Piezoelectric Disk 


Mode 

Natural Frequencies (cps) 

Experimental 

MARTSAM 

Mesh 

NASTRAN 

Mesh 

1 

22042 

• .23298 

24323 

2 


59805 

61114 

3 


103048 

104689 
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Figure 1 — Finite Element Mesh of Ferromagnetic Sphere; 
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IMPROVEMENTS IN SPARSE MATRIX OPERATIONS OP NASTRAN 


Shinichiro Harano 
Hitachi, Ltd. 

SUMMARY 


This paper describes improvements in sparse matrix operations for 
the NASTRAN program achieved by Hitachi, Ltd. (Japan). To solve a large 
scale problem at a high speed, a great emphasis was laid on how to make 
a reduction in execution time needed by matrix operations, since! the size 
of the problem depends largely on speed of matrix operations as well as 
on hardware and program performance. The descriptions in this paper are 
presented under Introduction plus five subjects: Sparse Matrix and Matrix 
Packing, Matrix Decomposition, Forward Elimination and Backward 
Substitution, Eigenvalue Extraction Methods and Parallel Processing 
Oriented Matrix Operations. These improvements can be applied to other 
versions of NASTRAN with a slight modification by using several 
subroutines which we have developed. 


INTRODUCTION 


Since the introduction of NASTRAN level 15.5.1 in 1974, we have 
improved it by a series of program enhancements. Highlights of them are 
development of the IG/OG (input Generator/Output Generator) program to 

perform automatic meshing and edit the results of calculation, a nd 

addition of isoparametric elements of two dimensions, three dimensions or 
axi-symmetry . Dealt with in this paper is another highlight of them. 

Recent drastic improvements in hardware performance have brought 
a gradual moving from the third generation computers, typified by IBM 570 
to distributed computers, also typified by IBM 5033. Being in step with 
such worldwide trends, Hitachi has developed HIT AC M-180 closely 
comparable with IBM 570/168 and HITAC M-200H providing a throughput three 
times that of IBM 570/l68. Along with these hardware breakthroughs, the 
parallel processing feature appearing with vector and array processors 
will be increasingly brought in, changing a current software environment 
greatly. 

Accordingly, in putting a further refined processing system flor 
matrix operations into practice under such situation, one must direct his 
attentions to hardware as well as software dimensions of the break- 
throughs. We have confirmed that a more effective use of a vector 
processor is well attainable by the Gaussian elimination of inner product 
type, also called "the Skyline method", which was proposed by Prof. 

Wilson (UCB:) , rather than by the conventional band matrix algorithm. 
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SPARSE MATRIX AND MATRIX PACKING 


Generally, matrices for structural analysis are characterized by 
sparsity. To take full advantage of this characteristic in matrix opera- 
tions, NASTRAN carries out matrix data packing. The way of matrix data 
packing is especially important for a problem where a large scale matrix 
is to be handled efficiently. So far, in transmission of matrix data 
between a secondary storage device and a main memory via an input/output 
buffer, packing routines have been used to transmit data from the input/ 
output buffer to allow matrix data to be referenced. However, this 

method is less advantageous to handle a large volume of matrix data since 

it needs much overhead time for the transmission. 

New "non-transmit" packing routine has been added to our version of 
NASTRAN to allow matrix data to be refered to directly from the input/ 
output buffer. The matrix packing format obtained as a result is shown 
in figure 1. In this figure, the string is a set of successive non-zero 
terms, plus the row number and length of string a hea d of t he se ter ms . 

Further, the number of strings in one column that are resident at one 

input/output buffer is given to control how to refer to matrix data resi- 
dent at the buffer. Padding information is also given to adjust a word 
boundary for data provided in double precision. At the present, the new 
packing routine to perform a direct reference to the input/output buffer 
makes only the READ option effective. This non-transmit type of routine 
called string by string receives or sends: 

(1) the start adress of a string in a buffer 

(2) the foremost row number of a string 

(3) the length of a string 

( 4 ) the instruction to show whether or not EOL (end of column detec- 
tion is to be made.) 

(5) type of matrix data 

Use of the packing routine permits various routines for matrix hand- 
ling to perform a direct reference to the input/output buffer if once 
they have received data addresses. The packing routine offers a buffer- 
by-buffer backspace feature for efficient backspacing in sequential 
access. Unlike a conventional backspacing that needs twice back record 
for a single read of one record (one column), as shown in figure 2, this 
feature omits overlapping of READ operation and back record, as also 
shown in figure 3* This feature eliminates the necessity of writing, in 
decomposition of a symmetric matrix, of a portion of the matrix to its 
upper triangular matrix from the last to the first columns of the sym- 
metric matrix, thus saving time for generating the upper triangular 
matrix. Furthermore, the feature requires the writing of only a lower 
triangular matrix onto the secondary storage device, bringing 10 to 30 % 
reduction in use of the disk space of the storage device. 

This new matrix packing technique is fully employed in the matrix 
decomposition described in MATRIX DECOMPOSITION. Figure 4 reveals how 
the technique is superior to conventional techniques by comparisons of 
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packing routines to pack and unpack one column of the matrix in respect 
to CPU time versus non-zero term densities. The CPU time for packing/ 
unpacking of one column is on the ordinate, while non-zero term densities 
are on the abscissa. The length of one column is 2000 and the whole CPU 
time has been obtained as the result of 300 iterations. Packing routines 
mutually compared in this figure ares 

(1) INTPK ' : Clears the area for one column to zero to perform ele- 

ment-by-element unpacking. (This is a conventional 
unpacking routine.) 

(2) UNPACK : Unpacks a column at one time. (This is a conventional 

unpacking routine.) 

( 3 ) PACK : Packs a column at one time. (This is a conventional 

packing routine.) 

( 4 ) INPNT : Clear the area for one column to zero and transmit 

string data directly from the buffer to the appropriate 
location of the column. (This is the new packing 
routine.) 

Figure 4 also shows CPU time for READ and WRITE operations in case 
of GINO (General Input Output Routines) as additional information. 
Although the READ and WRITE operations may be performed irrespectively of 
non-zero term density of the matrix, they cannot take advantage of spar- 
sity in case of low density due to unsatisfactory efficiency. Further, 
all elements including non-zero ones are written out onto the secondary 
storage device, making an increase in disk storage space needed. 

As a result, we adopted a combination of the clearing a core space 
to zero and the new routine of non-transmit type to unpack columns in 
matrix decomposition. In short, figure 4 reveals that the new unpacking 
routine permits a speedup approximately 2.1 times the conventional . 
unpacking routines in such a situation that densities of unpacked one 
column usually fall in the range of 0.4 to 1.0, owing to the method of 
matrix decomposition described' below. 


MATRIX DECOMPOSITION 


In solving a given system of linear equations, where the coefficient 
matrix is a large scale sparse matrix, it is general to decompose the 
matrix as, 

■ A = L*D*U (l) 

where, A is the original matrix (e.g. stiffness matrix), L is a unit 
lower matrix, U is a unit upper matrix and D is a diagonal matrix 
which is usually part of the diagonal portion of L or U . The discus- 
sion below assumes that the original matrix -is symmetric. To decompose 
the matrix, the Skyline method was used. The name of "Skyline" is 
derived from the fact that the contour line of column's all foremost 
none-zero elements is similar to a skyline. This method divides a por- 
tion enclosed in a skyline and a diagonal line into two groups whose 
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sizes are such that the groups are capable of being resident at a free 
area of a virtual djtorage space. The contents of these groups may be 
read from the secondary storage device which the results of calculation 
may be written out onto as needed. This is shown in figure 5. 

.1 Unlike conventional techniques, the Skyline method employs an upper 

i triangular matrix. The diagonal matrix is generated on a diagonal terms 
S iof triangular matrix U . The original symmetric, matrix is decomposed 
§ Jin the following algorithm. 




\ A— -f 
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The method carries out (2) through (4) in succession for j=2,...,n 
where, n is the dimension of matrix A . d^^ = a^^ is assumed. First, 

. the algorithm for the Skyline method used for an incore routine is desc- 
ribed with the help of the explanatory illustration of figure 6 . This 
method employs the Gaussian elimination of inner product type, as shown 
in this figure. All column's elements from foremost non-zero ones to 
f diagonal ones, including zero elements, are stored on the memory. The 
c, store address at that time is resident at array M. The address of the 
c ! I-th row diagonal term is expressed by M(l). The value at the I-th row 
0 and the J-th column (position pointed by Ij) is determined by the inner 
, product of vectors P and Q . Vector P , has a length of: 
i 

3 | JH = M(J) - M(J-1) (5) 

3 If JH=1, the diagonal terms are those obtained by matrix decomposition, 
and therefore, processing is skipped. Vector Q a party of inner pro- 
, :duct with P , satisfies: 

rt 

r J > I > J-JH ( 6 ) 


Let HT be a length of the intersection of vectors P and Q . Since 
the length of vector Q is, 

; IH = M(l) - M(l-l) (7) 

length NT is, ' 


NT = Min { I - (J - JH), IB’}L 1 


The following processing is performed only if NT is positive. Let 
NS be a start address of element being involved in inner product calcu- 
lation of vectors P and Q . Then, NS may be expressed as: 



NS = M(l) - NT 


(9) 


The address of the I-th row and the J-th column element to which the 
inner product value is to be added is, 

IJ = M(J) - (J-I) (10) 

The distance (IC) between the addresses of foremost elements of vectors 
P and Q is, 

IC S IJ - m(i) (11) 

Thus, letting X be a vector storing elements on the memory, the foremost 
elements for inner product calculation of vectors Q, and p are, . 

x(NS) and X(NS+IC), respectively. (12) 

The length of inner product, then, is NT, Assuming that the value of 
inner product between vectors P and Q is S, the I-th row and the 
J-th column is obtainable by: 

X(IJ) = X(IJ) - S ( 13) 

By performing the calculation for all I restricted by (6), the 
entire J-th column may be obtained. Notice that the value obtained by 
(13) corresponds to that by (2), In a practical LU-decomposition, as 
shown by (3), each element of the J-th column must be divided the corres- 
ponding diagonal,,! element (See (15 )). The J-th row and the J-th column 
diagonal term is, 

x(ij) . = x(ij) - ;f; ■ ' (l=J, kd=m(k)) (14) 

K=NS 

The last result for the I-th row and the J-th column is, 

X(K) = (K=NS , ..., NE) (15) 

By carrying out the above process for 2^J^n, the upper triangular matrix 
of matrix A may be generated in X, 

The above algorithm is well applicable to matrix calculation if all 
matrix elements are capable of being stored on the memory. Otherwise, 
matrix grouping is needed prior to implement the Skyline method. 

Assume that the original matrix is such a matrix as shown in figure 
7. First, this matrix is divided into some groups, each of which should 
not have more elements than those restricted by a core space. In the 
example of figure 7» the matrix is divided into four groups, each of 
which has not more than 15 elements. These divided groups provide the 


18 



information of table I f headings of the table are: 

(1) K(l) : the~ current group number 

(2) K(2) : minimum group’ number needed by calculation of group K(l) 

(3) K(3J s minimum column number in group K(l) 

( 4 ) K 44 ) 5 maximum column number in group K(l) 

(5) M(J) : pointer array of the diagonal term in the J-th column 

With the matrix grouped above, the algorithm of the Skyline is pro- 
ceeded as follows. Symbols used in the description ares 

X : open core array 
I A : start address of group A element 

IB : start address of group B element 

IM : start address of pointer array of diagonal terms 
KA(l), KB(l) : group numbers for groups A and B 
KA(2), KB(2; : minimum group numbers involved in calculation of 
groups A and B 

ka(3), kb ( 3 ) s minimum column numbers of groups A and B 
KA(4), KB(4) s maximum column numbers of groups A and B 
Apart from the explanation of how to determine these values, we begin 
our discussion with the following assumptions. The areas are already 
assured for groups A and B (headings are X(lA) and X(lB)) and for the 
address of diagonal terms (heading is X(lM)). The diagonal term address- 
es are set in advance. Let NEQ be the number of unknowns of a given 

system of linear equations, and NGP be the number of groups. All control 

information for NGP groups, except those for diagonal term addresses, are 
generated on a scratch file in advance. 

Then, the following steps are repeated for P until NGP. 

P = 1, ... , NGP (P J group number in current calculation) 

/ 

Mutually permutated, for Pj*l, are the address of A and that of B, and 
group information of A and that of B, that is, 

IAS^IB, KA(k)i^KB(k) (k=l,...,4) (16) 

Let Q be the minimum group number needed to generate group P, Then, 

Q = KA(2) 

For PjisQ , group Q is already on B if Q = P-1; otherwise the control infor- 
mation about group Q is read from the scratch file to be set to KB(k) , 
where ,k=l , . i ,4 . 1 Then, the correlation values between groups P and Q are 
added to group P. This process is carried out for Q by an increment of 1 
until Q = P-1. For Q = P, the correlation between P and Q becomes that 
between P and itself on A. After the completion of the above step for 
group P, the elements of it are written out onto the secondary storage 
device. The implementation of these procedures for individual groups 
(l-P-NGP) allows an upper triangular matrix to be generated in the column 
direction on the secondary storage device. 

The correlation between groups P and Q is obtained as follows. 
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Assume that group P is on A (heading : X(lA) area), while group Q is on 
B (heading : X(IB) area). Let NTA be the number of columns of group P 
on A, and NTB be the number of columns of group Q on B. Then, 


NTA = KA(4) - KA(3) + 1 
NTB = KB ( 4 ) - KB(3) + 1 


(17) 

(18) 


Similar to figure 6, handling of groups P and Q is shown in figure 8, 
Let J=l, ... , NTA be column numbers of group P on A, Then, the column 
number of P on the entire matrix is, 


JJ = KA(3) + J - 1 
The length of column J is, 
JH = M(JJ) - M(JJ-1) 
where, if JJ=1, 

JH = M(JJ) 


(19) 

( 20 ) 
( 21 ) 


Similarly, let 1=1, ... , NTB be column numbers of group Q on B. Then, 
the column number of Q on the entire matrix is. 


II = KB(3) +1-1 
The length of column I is, 
IH = M(ll) - M(II-1) 
where, if 11=1, 

IH = M(ll) 


( 22 ) 

(23) 

(24) 


Again, let NT be the length of inner product of column J in group P and 
column I in group Q. Then, 


NT 


= Min | IH, JH-(JJ-Il)j -1 


(25) 


Also, let NS be the start address involved in the inner product of 
column I in group Q and NE be the address of the last portion. Then, 


NS = M(II) - NT 
NE = M(II) - 1 


(26) 

(27) 


The displacement (IJ), where the inner product value is added on area A, 
is expressed as: 


IJ = m(jj) - JJ + II 


(28) 
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If NS>NE, column I in group Q, has nothing to do with the calculation for 
' column J in group P; otherwise, inner product must be calculated. Since 
v the start addresses of areas A and B on the open core are I A and IB, 
respectively, by putting, 

o IC = IJ - M(II) (29) 

( 

8 ithe inner product is, 

q •; 

i~ NE 

3. - X(IJ) = X(IJ) - Y_ X(IB+K-1)*X(IA+K+IC-1) ; (30) 

12 K=N8 

' The calculation of contributions from group Q to group P is completed if 1 
3- the above steps are carried out for all columns in group Q. 

;.c 1 . 

1 ? ‘ After calculations on all groups such that, 

ka( 2) — Q — P—1 (31) 

'?.■ the autocorrelation of group P itself is obtained by the same procedure 
22 as that used by previous calculation of the correlation between groups P 

2 ./ and Q by regarding that area A is the same as area B. In this calcula- 

tion, I varies within the ranges 

'y-\ ; 

y • 

2oj 1<I<J (32) 

If I=J, (30) must be replaced by the procedure below. For K such that, 

■ r> ! 

| 

.30 ! NS<K<NE (33) 


>2 ID = IJ - KJ + 1 is obtained. If KD = M(lD) is established, the column J 
X> ;is completed by (14) and (15), 


33 The Implementation of the above procedure for the range of (32) 

gives the autocorrelation of group P. The results are written out onto 
> an upper triangular matrix data-block for each column. 


The interchange of addresses and control information by (l6) are 
•j needed for handling a succeeding group. This takes the place of moving 
■- 1 group P completed on area A to area B simply by interchanging addresses 
: r ' and control information. This eliminates a data transmission, and fur- 
-'3 jther allows one group to be read in the core only once if the correlation 
22 between two groups affects only adjacent groups. 

: 


2 As already suggested, the matrix must be prepared for grouping prior 

to the implementation of the algorithm of the Skyline method if matrix 
data overflows the core space available. First, the size of each group 
r must be determined to provide such control information as in table I , 

The size of a group depends upon how the open core is large at execution. 
Figure 9 presents an open core layout. Given the open core size (NX) , 
the size of the memory to be allocated to one group may be obtained by 
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bisectioning the axea excluding a working area needed. This allows up 
to MAXT elements of one group to be stored. By use of the new unpacking 
routine, only start row numbers of each matrix column are picked up to 
create table I information. Such information are stored on a scratch 
file in such a way that each group is on one record. At the same time, 
the addresses of diagonal terms of each group at a working area are also 
stored on array M. Thus, the preparation is completed. 

To read group P, The unpacking method must be used in which the 
input/output buffer can be directly re fer red to from an input matrix 
data-block. This is also applied to reading Q, where the buffer is 
directly referred to from an upper triangular matrix data-block. After 
the completion of calculation, each group is written out following the 
end of the upper triangular matrix. The diagonal elements are also 
written out onto another output data-block for succeeding forward elimi- 
nation and backward substitution. 

The following are the results obtained by applying the Skyline 
method to practical examples. Table H gives matrix characteristics for 
four data with a comparision of matrix characteristics in case of the 
conventional band matrix method. Figure 10 compares CPU time for the 
band matrix method with that for the Skyline method. The Skyline method 
in these examples gives two cases in which a vector processor has been 
applied and it has not been applied. The vector processor has been also 
applied to the band matrix method, resulting in no improvement of CPU 
time. This figure, therefore, does not that case. Figure 10 reveals 
that the Skyline method consumes 33 to 66 % of CPU time needed by the 
band matrix method. If the vector processor is applied to the Skyline 
method, this value drops to as many as 16 to 28% of the CPU time. Fur- 
ther, for data of 7000 degrees of freedom, the Skyline method has needed 
CPU time on the same percent basis as the band matrix method. 


FORWARD ELIMINATION AND BACKWARD SUBSTITUTION 


NASTRAN is designed to proceed forward elimination and backward 
substitution while retaining vectors of load terms on the memory as much 
as possible. As the result of matrix decomposition by the Skyline 
method, an upper triangular matrix is generated and diagonal terms are 
stored on another file. This means that forward elimination is an inner 
product type and that store-type calculation is to be carried out after 
division of load terms by diagonal elements. The procedure for solving 
a system of linear equations: 

uSuX = B (34) 

is separated into the forward elimination process 

U^Y = B (35) 


and the backward substitution process 
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FX = D _1 Y 


(36) 


The forward elimination process is shown in figure 11. In this figure, 
1. . ? forms a string starting with the i— th column and the 

oi’lU upper triangular matrix (string length is 3 here). 

f For this string, the following calculation is carried out. 

^ r 



9 


b. = b. 
1 a 1 a 


j+ST-1 

-X ilk * *ka 


(a«l,2,3) 


(37) 


; This is performed for all i-th column strings in the upper triangular 
matrix. Then, this procedure is repeated after setting i=i+l, Since 
, '(37) is an inner product type, high speed calculation is possible. 

Further, as for string’s elements, the new packing routine refers to 
_."t the input/output buffer,, . saving time for data transmission and reducing 
'■/;< time for calling the unpacking routine due to string— by-string call un- 
like conventional element-by-element call. Thus, the forward elimination 
: :0 process is inner-product-type operation for matrix's factors decomposed 
; by the Skyline method, while the forward elimination process of the con- 
ventional method is store-type calculation. 

On the other hajjid, the backward substitution process begins with 
the generation of D” X Y for load term Y already generated by the forward 
2C 'elimination process. The process allows calculations independent of the 
following process because all diagonal term elements are already genera- 
A ted on a file as one vector on matrix decomposition. After the division 
y; of load terms by diagonal terms, backward substitution is performed by 
y: (backspacing the upper triangular matrix file buffer-by-buffer. This is 
shown in figure 12. In the process, 

2. y ka = y ka “ u kj* x ja (k=i,i+ST-l; a=l, 2, 3) ( 38 ) 

is calculated for each string of each upper triangular matrix column. 
y- At that time, the non-transmit unpacking routine is called for each 
•y- string, and only addresses in the input/output buffer are passed to the 
;*■ routine of forward elimination and backward substitution. 

- . Figure 13 gives the results of forward elimination and backward sub- 
stitution by use of data shown in table II :. This figure reveals that | 
y the new method requires osly.16 to 54% of CPF time needed by forward eli- 
mination and backward substitution of the conventional method. Taking 
Hi, it into account that the coding for both processes are not oriented to 
; - the vector processor, more satisfactory results will be expected in res- 
■f. ipect to CPF time by further improvements. 


EIGENVALFE EXTRACTION METHODS 


In real eigenvalue extraction methods we attemped to develop these 
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methods in two directions: that is, partly the speeding up of the Inver- 
se Power method, partly the development of a simultaneous iteration me- 
thod. — 

At first we describe the Inverse Power method. Eigenvalue extrac- 
tion methods are generally divided into two groups: tracking methods 
( the Inverse Power method and the Determinant method ) and transforma- 
tion methods ( the Householder method and the Givens method ) . Though 
transformation methods are able to solve rapidly an eigenvalue problem 
in the range of comparatively small scale problems, large scale problems 
are unfavourable to them. On the contrary, the Inverse Power method has 
been employed frequently owing to less restriction than transformation 
methods . 

Since the Inverse Power method in NASTRAN is accompanied by move- 
ments of shift points, it needs to use iteratively matrix decomposition 
and PBS (forward elimination and backward substitution). Improvements 
that were previously mentioned were applied to the Inverse Power method, 
so that we could Improved the CPU performance of the Inverse Power me- 
thod which is two or three times as much efficient as that of the con- 
ventional Inverse Power method. Table HI shows that the CPB performance 
of the new method without the vector processor amounts to 2.7 times that' 
of the old one. If the vector processor is applied to the former, the 
CPU performance of it will be equal to about 3.5 times that of the con- 
ventional one. 

Now we describe a simultaneous itetation method which is called the 
Jennings method. There is a problem to find q eigenvalues in ascending 
order from the lowest value and q eigenvectors corresponding to them 
for the general eigenvalue problem: 

K x = r M x (39) 

where K is a symmetric matrix of positive definite type and M is a 
symmetric matrix of non-negative definite type. The Jennings method is 
useful for calculating a set of eigenvalues from the lowest value and 
has no weakpoint that some important eigenvalues are often missed in 
calculation. The algorithm is shown in figure 14. This method is diffe- 
rent from the Subspace Iteration method on operations of orthogonaliza- 
tion shown in (i), (j), (k) , (l). 

The Jennings method needs the following input data: 

(1) ND : number of eigenvalues to be extracted (2 — ND290). 

(2) LMAX: the maximum niimber of subspace iterations (the default 

value is l6)Y 

( 3 ) IEP : convergence parameter (if IEP<0, then EPS=10 ; 

otherwise EPS=0.000l) . 

The dimension "m" of subspace is decided on by 

m = min(n, 2q, q+8) (40) 

The selection of initial iteration vectors is most important for the 



convergence of subspace and the convergence ratio is decided on by the 
"neighborhood" between subspace spanned by eigenvectors and subspace 
spanned by initial iteration vectors. Assume that 


K = (kj^ , K = (-„) 


: k.. =k.. , m. . = m . . 
ij Ji 13 ji 


(41) 


Then, k^ £ 0 is always satisfied as K is positive definite. The 
matrix Go which is composed of m initial vectors will be generated 
as follows; 

(1) At first, the first column of Go is a vector D, where 
D(I)= m. . 

(2) check k.,^0 » and 


D(I) = D(I) / k u = m u / k u 


(42) 


<3) 


sort D(l) (1=1, ... , n) and select (m-l) values in decending 
order from the largest: 


d(i 2 )2d(I 3 )2...£d(I u ) (43) 

where the i-th vector of Go (i= 1, ... , m) is the unit vector, the i-th 
component of which is equal to 1. 


The criterions of convergency are as follows: 

(1) q eigenvalues and q eigenvectors are extracted. 

(2) number of subspace iterations amounts to LMAX. 

( 3 ) there is no CPU time to execute three subspace iterations, 
because output of the results needs a little time. 

Let the eigenvalues in the L-th iteration loop be on a vector E(l) (l= 1, 
... , m) . After reordering E(l) in ascending order, the eigenvalues in 
the (L-l)-th iteration loop are stored on a vector E*(l) (l= 1, ... , m) , 
If E(l) satisfies the following relation* 


E(l) - E*(l) 
E(I ) 


< EPS 


(44) 


then E(l) is already converged; otherwise E(l) is not converged. If (44) 
holds true for all I (lSI<q), the convergence will be achieved owing to 
criterion (l) . If (44) doesn|t hold true for some I and number of sub- 
space iterations is equal to LMAX, the calculation of eigenvalues will be 
stopped due to criterion (2) . 


The orthogonalization of subspace is also important for a simultane- 
ous iteration method. If a mass matrix M is non-positive definite and 
operations of orthogonalization isn't applied to subspace during itera- 
tion loops, the orthogonality of subspace will be breaking. For a eigen- 
value problem with a non-positive definite mass matrix, the orthogonali- 
zation of subspace is necessary for iteration vectors to converge to 
eigenvectors. Consequently we adopted the Jennings method which ortho- 
gonalizes iteration vectors just after the calculation of eigenvalues in 
subspace. The generalized Jacobian method is adopted in the eigenvalue 
extraction on subspace. 


25 



Tablell shows that the CPU performance of the Jennings method with- 
out a vector processor (or with it) is 4*0 (or 4.l) times as high as that 
of the conventional Inverse Power method. Thus, the Jennings method 
consumes about two-third of CPU time used by the new Inverse Power meth- 
od, This fact results from the following reason. While the Inverse 
Power method needs several decompositions of full size matrices in every 
movements of shift points, the Jennings method is more efficient to solve 
.large scale eigenvalue problems than! the Inverse Power method, for the 
former needs only one decomposition of full size matrix and several deco- 
mpositions of small scale matrices oh subspace. 


PARALLEL PROCESSING ORIENTED MATRIX OPERATIONS 


Our vector processor adopts a pipeline system and uses a compiler 
system in which a FORTRAN source program is translated into a set of 
instructions specially for the vector processor with recource to the 
option active in compilation. This means that the object program gene- 
rated by the compiler depends on the skillfulness of coding. 

To disouss more specifically, this section presents the results, of 
our test. In this test, we measured CPU time per single term for the 
length of a DO loop (string length) by carrying out thiree inner product 
type operations and one store type operation, in order to determine the 
basic operation in matrix decomposition. The results are shown In figure 
15, As for the inner product type operations, two cases were further 
considered: the case where the vector processor was applied and the case 
where it was not applied. The examples of coding 1 ! used in our test are: 


(l) Complete inner product type 


REAL* 8 A(1000), B(l000) , X(lTER), SS 

DO 10 I = 1 , ITER 
SS = 0.0D0 

DO 20 J = 1,LL 1 

20 SS = SS + A(J)*B(J) 1 

X(l) = X(l) - SS 
10 CONTINUE 


Inner product loop 


(2) Index explicit type 

REAL* 8 X(l), SS 

DO 10 I = 1, ITER 
SS = 0.0D0 

DO 20 J = 1,LL j Explicit index type 

20 SS = SS + X(lA+J_i)*x(iB+J-l) 1 inner product loop 

X(IC+I-1) = X(IC+I-1) - SS 
10 CONTINUE 
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(3) Subroutine inner product type 

i 

REAL* 8 X(;l) , SS 

DO 10 I = 1, ITER 
SS = 0.0D0 

5 CALL DOTP (X(ia), ; x(lB), SS, LL) <■ 

X(IC+I-1) = X(IC+I_1) _ SS 
& _ 10 CONTINUE 

o ? 

s : 


(4) Store type (the three terms operation) 

REAL*8 X(l) ,SS 
DO 10 I = 1,ITER 
DO 20 J = 1,LL 

X(IC+J-1) = X(lC+J-l) + X(lA+J-l)*X(lB+J-l) 

>,; • 20 CONTINUE 

2 r 10 CONTINUE 

;x> The above examples of coding are only for our test and there is no mean- 
Vi ! ing in operation itself. The index ITER is the number of iteration loops 
and ITER = 2000 in our test. Though the store type (the three terms) 
operation is applicable to the vector processor by changing its indices, 
; \;at that time we left it as it was, and then it was not applicable to the 
vector processor. 

>3 

A close observation of figure 15 first exhibits that, as for the 
complete inner product type operation, use of the vector processor brings 
y; about an improvement in 'speed as much as 5.5 times that obtained by the 
same type of operation without the vector processor. Unfortunately, 
however, NASTRAN is not oriented to the way of coding for the complete 
inner product type, since it uses an open core as a working area. Thus, 

-s two possible ways for coding are explicit index and subroutine inner 
y: product types. Without a vector processor, the subroutine inner product 
LL 'type is more advantageous than the explicit index type. On the other 
Ly .hand, with the processor, the former is less advantageous than the y 

i y latter. 

Further, figure 15 reveals that, with the vector processor, the 
explicit index type almost keeps in step with the complete inner product 
type in respect to CPU time. However, without the processor, the former 
has consumed CPU time as much as 2,2 times that the latter has consumed. 
For a longer inner product loop, the subroutine inner product type 


SUBROUTINE DOTP: (A, B, SS, LL) 
REAL* 8 A(LL), B(LL), SS, S 
S = 0.0D0 
DO 10 I =s l.LL 
S = S + A(I)*B(I) 

10 CONTINUE 
SS = S 
RETURN 
END 
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without the vector processor is more advantageous than the explicit index 
type without it. This is attributable to that there is a difference in 
optimization level between the operation with the vector processor and 
that without it. The subroutine inner product type is advantageous if 
the subroutine's overhead time can be overridden due to a long DO loop; 
otherwise, it is less advantageous than other types in speed. The result 
of the store type operation is also exhibited in this figure; this type 
does not enjoy the maximum benefits of optimization. 

An observation of figure 15 also shows that the extent of optimiza- 
tion in various types of operations depends largely on program coding. 

Of course, it is ideal that the maximum optimization is always possible 
for any type of operation; however, the extent of optimization varies 
depending upon type of FORTRAN, Accordingly, in coding the algorithm for 
the Skyline method of matrix decomposition, we adopted the explicit index 
type if use of a vector processor was possible; otherwise, we used the 
subroutine inner product type. In the future we intend to use the expli- 
cit index type as long as the optimization feature of FORTRAN is satis- 
factorily refined. 

So far, the inner product type has been more advantageous than the 
store type in respect to speed thanks to use of registers. However, the 
advent of a vector or array processor is changing this situation. 
Actually, in case of HITAC M-200H, the latter has displayed almost the 
same performance as the former. Further improvements of the parallel 
processing systems may reverse the superiority of the inner product type 
to the store type. 

As already described, CPU time needed for the store type, inner 
product type operations accompanied with or without data transmission 
depends largely on how to make a program. Use of a higher speed computer 
and parallel processing system is greatly expected to change a current 
software environment to a large extent. Technological breakthroughs of 
software and of hardware would interact more closely in improving sparse 
matrix operations. 


CONCLUDING REMARKS 


In this paper we discussed about improvements in sparse matrix 
operations of NASTRAN. Recent advance of parallel processing systems has 
been changing surroundings in software. Especially, a vector processor 
attached to a general-purpose computer is favorable to a long DO loop 
operation. For example, the Skyline method which we have developed this 
time in the field of matrix triangular decomposition conforms to the 
pipeline control feature observed in the vector processor. On the cont- 
rary, the conventional band matrix method or the wavefront method which, 

adopt store type operations don't adapt themselves to the pipeline^. j; 

control system, for they need complicated indices operations and are 
difficult to deal with a set of arithmetic data as vectors. 
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What is more, the way of packing/unpacking and the method of forward 
elimination and backward substitution were conformed themselves to the 
Skyline method, so that the CPU time for solving a problem was reduced by 
half. Further, in real eigenvalue extraction we have improved the CPU 
; performance of the Inverse Power method and added the Jennings method to 
~ NASTRAN. The Jennings method is more effective in many cases than the 
' new Inverse Power method. I 
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TABLE I.- GROUP CONTROL INFORMATION 


Group 

K(l) 

K(2) 

K(3) 

K(4) 

M(J) 

Group 1 

1 

1 

1 

6 

B 

3 

6 

8 

12 

15 

n-r»r\n -r^ 0 

y x 4— 

2 

1 

7 

9 

B 

7 

14 



' 

Group 3 

3 

2 

10 

11 

B 

8 

B 


■ 

B 

Group 4 

4 

1 

12 

12 

li 




■ 

B 


K(l) : Current group number 

K(2) s Minimum group number needed by calculation 
of group K(l) 

» K(3) s Minimum column number in group K(l) 

K(4) s Maximum column number in group K(l) 

M(j) : Pointer array of the diagonal term in the 
J-th column 

, < ! 1 ~ r : ^ y *• jti. p . ■' 
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38 
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1024 
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477 
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TABLE III. CPU TIME'- OLD, NEW INVERSE POWER METHODS 


AND JENNINGS METHOD 


Data name 

Total degree 

of freedom 

• 

Extraction method 

• 

Vector processor 

Number of eigenvalues 
to be extracted 

CPU time (sec) 

Total 

Eigenvalue 

extraction 

Data 1 

208 2 

Old Inverse Power method 

Can not be applied 

1 0 

1864 

1813 

Data 1 

2082 

New Inverse Power method 

Yes 

1 0 

344 

293 

Data 1 

2082 

Jennings method 

Yes 

1 0 

241 

■I 

Data 2 

mm 

Old Inverse Power method 

Can not be applied 

1 0 

1 1 05 

1062 

Data 2 

2 1 27 

New Inverse Power method 

No 

1 0 

398 

345 

Data 2 

BB 

New Inverse Power method 

Yes 

1 0 

332 

278 

Data 2 

2127 

Jennings method 

No 

1 0 

276 

232 

Data 2 

2127 

Jennings method 

Yes 

1 0 

266 

222 
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FIGURE 5. 
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FIGURE 8. SKYLINE METHOD WITH GSiUPING 
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FIGURE 10. CPU TIME : SKYLINE METHOD AND BAND MATRIX 

METHOD 
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FIGURE 1 4 . ALGORITHM FOR JENNINGS METHOD 
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SOLUTION SENSITIVITY AND ACCURACY STUDY OF 
NASTRAN FOR LARGE DYNAMIC PROBLEMS 
INVOLVING STRUCTURAL DAMPING 

A. J. KALI NOWS K I 
NAVAL UNDERWATER SYSTEMS CENTER 

SUMMARY 

This paper is concerned with both the solution sensitivity and solution 
accuracy of large dynamic problems involving NASTRAN SOLUTION 8 (i 0 e„, the 
steady state dynamic response option wherein all response quantities vary as 
e iwt } where o> is the driving frequency and t is time). Using a submerged 
steel plate with a viscoelastic layer as the bench mark sample, the solu- 
tion sensitivity and solution accuracy is checked. The solution sensitivity 
is examined by running the same finite element model on different computers, 
different versions of NASTRAN, and different precision levels. The solution 
accuracy is evaluated for these same runs by comparing the NASTRAN results 
with the exact solution of the same problem. 

SYMBOLS 

[B] Damping Matrix 

Ci Dilational Wave Speed in Fluid 

{F} Applied Force Vector 

[K] Stiffness Matrix 

[T] Modified Complex Stiffness Matrix 

[M] Mass Matrix 

k Wave Number (w/ci) 

P.j Incident Fluid Pressure 

P Q Plane Wave Amplitude 

Pg Back Side Fluid Pressure 

P<j Front Side Fluid Pressure (Scattered Component) 

t Time 

{U} Solution Displacement Vector 
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SYMBOLS (Cont'd) 


x 


U) 



r 

»y 

i 

>p 


p 


Spatial Coordinate 
Driving Frequency 
Residual Solution Vector 
Real Elastic Lame 1 Constants 
Corresponding Viscoelastic Constants 
Material Mass Density 


INTRODUCTION 


This paper is concerned with the solution accuracy of 1, 2, or 3-dimen- 
sional steady state (time harmonic) structural and/or continuum problems whose 
response quantities all vary in time in proportion to e lu)t . The linear equa- 
tions of motion for such problems usually reduce to an expression of the form 

[-o) 2 [M] + io)[B] + [K] ] {U} = { F} (1) 

[K] 

where [M], [B], and [K] denote the mass, damping and stiffness matrices (MDD, 
BDD and KDD using usual NASTRAN DMAP notation), u> is the driving frequency and 
{ F> are the applied forces. The results presented in this paper focus on con- 
tinuum type (e.g., figure 1) applications with structural damping, however, 
once the form of equation (1) has been constructed, the solution becomes a 
matter of solving large banded symmetrical systems of complex linear algebraic 
simultaneous equations. Clearly, such equations can also be the end point 
resulting from many other NASTRAN steady state formulations, either from direct 
structural formulations or from related fields through analogies. Thus compari 
sons of solution accuracy, run time, etc. can be viewed and interpreted in a 
more general vein than simply applying only to problems of the type depicted in 
figure 1. 

The motivation for this comparative study resulted as a consequence of 
obtaining some unexpected results on some solution 8 (steady state time har- 
monic rigid format) problems similar to the one shown in figure 2, except for 
the fact that the initial model had inclusions throughout the rubber thus mak- 
ing analytical solutions to the problem unwieldy. 



PARAMETRIC STUDY MODEL 


In order to better understand the accuracy limitations of the results of 
the initially more complicated inclusion filled model, a simpler homogeneous 
layered model (figure 2) was constructed and physically corresponds to a 
totally submerged 2.0" steel plate with a 3.05" viscoelastic rubber layer 
glued to the steel surface. The input corresponds to an incident pressure 
wave 

p i = p 0 el(kX + Wt) > k = w/c i 

where x is the horizontal coordinate along the line of propagation, c x is 
the dilatational wave speed in the fluid, and P 0 is the plane wave amplitude. 
The exact analytical solution to this problem is known (ref. 1), consequently 
an accuracy check on the finite element solutions is available. Clearly, the 
figure 1 model is a spatially one-dimensional problem, consequently the cor- 
responding finite element model need only be one element wide as was done, 
for example in ref. 1. However, the finite model was made up to eight ele- 
ments wide for the following reasons: (1) the model simulates the more com- 

plex model except for the fact that the inclusions are removed by filling 
their space are with uniform elements having the same material property as 
the surrounding rubber material; (2) the problem is artificially made mathe- 
matically larger so that more meaningful comparisons of CPU run times could 
be made; (3) larger problem sizes tend to draw out any potential problems 
with equation solvers. It is not our intent to discuss or explain the setup 
of wave propagation problems of the type represented by the figures 1-2 
example model; the reader is referred to refs. 1 and 2 for supplementary 
details. In fact, the demonstration problem used here is very similar to 
the one used in the ref. 1 sample problem except that the plate and visco- 
elastic thicknesses are different, the damping coefficient in the visco- 
elastic layer is different and that the steel plate is represented here 
approximately with CBAR elements rather than with solid elements as in ref. 1 
Specifically, the material constants employed are listed below 

DEMONSTRATION PROBLEM PHYSICAL CONSTANTS 


MATERIAL 

A*" 

psi 

r 

VI 

psi 

H 

i 

V 

psi 

1 b-sec 2 /in l( 

Water 

345;600. 

0.0 

0.0 

0.0 

.000096 

Steel 

17,307,000. 

11,538,000. 

0.0 

0.0 

.000735 

Viscoelastic 

Material 

86,703. 

115.9 

8670.3 

11.59 

.0003599 


where the meaning of the elastic and viscoelastic constants are defined in 
detail in ref. 1. 
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Since the topic of interest here is related to the class of problem 
treated by ref. 1, it appears appropriate to print an errata to the ref. 1 
paper 

• in equation (15) of ref. 1, replace 622 = x r with G22 = \ r + 2y r 

• in equation (16) of ref. 1, replace G22 = e with G22 = e(l + 2yV^) 

• in equation (2) of ref. 1, replace w 2 with -u> 2 

• in equation (17) of ref. 1, replace +n d with +in. in the k 2 defi- 
nition 2 2 


PARAMETER VARIATIONS 


The basic finite element model, figure 2, was exercised for a frequency 
sweep of 7 different incident frequencies (3.0 kHz, 4.0 kHz, 6.0 kHz, 8.0 
kHz, 17.5 kHz, 22.5 kHz, 35.0 kHz). Running the figure 1 model on NASTRAN 
for the above frequency sweep is designated as a typical run and correspond- 
ingly assign it a "run number", which runs from the number 2 through 9. Run 
number 1 is the exact solution and therefore is the only non NASTRAN designa- 
tion.^ is called a run since even the analytical solution involves a computer evalua- 
tion). 

Next, the same frequency sweep input data was rerun while varying the 
following parameters: 

• solution precision (S.P. or D.P. on the same computer) 

• type of computer (UNIVAC 1108; DEC-VAX; CDC Cyber 175) 

• level of NASTRAN (both NASA and MSC versions are considered) 

• date (i.e., the same input is resubmitted on the same computer, 
using the same version of NASTRAN but on different days) 

The last parameter (i.e., the date) seems a waste of computer time, however 
as is shown later, some unexpected results are encountered. 


DMAP INSTRUCTIONS FOR PRINTING SOLUTION ERROR RESIDUAL 


It is of interest to know the accuracy of the solution solving capa- 
bility of the equation solver used by the particular version of NASTRAN 
employed by the user. Specifically, if the solution (U) is found by NASTRAN, 
how well does it satisfy the linear simultaneous equations (1)? Consider 
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substituting the solution {U} into equation (1) and then transposing the 
applied force vector to the left hand side of equation (1) to obtain 

[K]{U> -(F) = (A) (3) 

If the equations have been solved exactly, then the residual vector {a} will 
be identically zero. The appearance of large nonzero entries in {a} would 
imply potential inaccuracies in the solution vector {U}. The question of 
"how large is large?" should be viewed by comparing the size of a particular 
entry in the {A} vector to the size of the applied loads (for this reason, 
the load vector { F} is also printed). For example, a residual of .2 would 
be a big residual if the applied forces are on the order of 1.0 lbs.; however, 
if applied forces are say 100,000 lbs., the .2 residual is acceptable. 


In order to print out the residual vector {A} for the 3,000 Hz driving 
frequency case, the following DMAP instructions were used (note, the {a} 
vector is printed with the heading DELS0L). 

• For UNIVAC 1108, SOL 8, LEVEL 17.0 (Runs 3a, 4a) 


ALTER 159 

ADD5 KDD,BDD,MDD,,/KBARX/C,Y,ALPHA=(1.0,0.0)/C,Y,BETA=(0. 0,18849. 5592)/ 

C,Y,GAMA=( -355305758. 44, 0.0) $ 

MPYAD KBARX,UDVF,PDF/DELS0L/C,N,O/C,N,1/C,N,-1/ $ 

MATPRN DELS0L,PDF, ,,//'$ 

ENDALTER 

CEND 

where frequency dependent input constants BETA and GAMA are to be i nput by 
the user and are simply defined as: 

BETA = 0.0 + io) 

GAMA = -a) 2 + i 0.0 

where u> = the driving frequency in radians/sec 

• For UNIVAC 1108, SOL 8, NASA LEVEL 15.5 (Run 2) replace ALTER 159 
with ALTER 139 

• For VAX, SOL 8, NASA LEVEL 17.5 (Run 5b) same as 1108, NASA LEVEL 
17.0 

• For CDC CYBER, SOL 8, MSC LEVEL 48B (Run 6) replace ALTER 159 with 
ALTER 139 

• For VAX, SOL 8, MSC LEVEL 52 (Run 7) replace ALTER 159 with ALTER 
139 

• For VAX, SOL 8, MSC LEVEL 60 (Run 8) replace ALTER 159 with ALTER 
139 

• For VAX, SOL 26, MSC LEVEL 60 (Runs 9b, 9c, 9d) 
replace ALTER 159 with ALTER 409 

replace UDVF with UHV (3rd line) 
replace PDF with PD (3rd and 4th line) 
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When a large residual is encountered, it is desirable to know the node 
and component number where a large residual appears (i.e., knowing the run 
number of the questionable residual, what node-component number does this 
correspond to?). By inserting a DIAG 22 card, the desired correspondence 
between run number of {A} and the node-component number can be made. 

It is important to note that the simple DMAP sequence as presented will 
apply to only a single frequency; thus, if a frequency sweep is employed, 
only the ntrr column of the DELSOL vector (i.e., {A}) will be correct; the 
remaining columns of DELSOL should be ignored, where n = the n^* 1 value in the 
frequency sweep appearing .on the MA5TRAN FREQ card. -ard. 


DISCUSSION OF RESULTS 


The primary variables of interest to us in this study are: (1) the 

transmitted pressure in the fluid on the back side of the steel plate, Pg, 
(e.g., in element number 100352 as shown in figure 2) and (2) the scattered 
pressure in the fluid on the front side of the plate, P$ (e.g., in element 
number 100378 as shown in figure 2). The transmitted pressure is read 
directly from the NASTRAN printout, whereas the scattered pressure is obtained 
indirectly from the NASTRAN printout by simply subtracting the incident 
pressure (equation (2)) from the total pressure printed by NASTRAN. The 
scattered pressure is of prime importance with regard to establishing the 
energy absorbing properties of the viscoelastic configuration. As discussed 
in ref. 2, it has been our experience that for steady state wave propagation 
problems of the type considered here, at least 10 elements per wave length 
are needed to adequately compute the pressure response for elements of the 
type employed in this study. In order to demonstrate this accuracy limita- 
tion, the model has been purposely exercised in a driving frequency range 
that is too high for the mesh to properly produce sufficiently accurate 
results (i.e., the mesh is too coarse for some of the higher frequencies). 

For a rubber wave speed of C 2 = 15540. in/sec, and the coarsest element in 
the rubber mesh (.1" x .1" elements), the following frequency vs. element/ 
(wave length) chart is constructed: 
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freq. (kHz) 

elements/wave length 
(rubber in figure 2 ) 

3.0 

51.8 

4.0 

38.8 

6.0 

25.9 

8.0 

19.4 

17.5 

8.8 

22.5 

6.9 

35.0 

4.4 


Based on the above chart, it is expected that the accuracy of the finite 
element solution in relation to the exact solution should start to drift at 
frequencies of 17 » 5 kHz and higher. 

A total of 14 computer runs were made (designated as runs number 2, 3a, 
3b, o o . o 9d) and are tabulated in Table 1 (for the transmitted back side 
pressure P B ) and in Table 2 (for the scattered front side pressure, P 5 ). 

The pressure results are normalized by the magnitude of the incident pressure 
wave, P 0 „ In addition to the pressure result. Table 1 has two additional 
pieces of information, namely the magnitude of the largest complex residual, 

I A 1 , appearing in the residual vector {A} as computed by equation (3), (the 
j A | value is obtained by scanning the DMAP printed DELS0L printout and seek- 
ing out the largest absolute value of all the rows of the {a} vector corres 
sponding to the 3,000 Hz frequency case. Note that the residual for only 
the 3,000 Hz case is reported. Also listed is the total CPU time required 
to execute the full frequency sweep solution for the run in question. 


Word Length (Precision) Sensitivity 


The earlier Level 15.5 version of NASA NASTRAN for a UNIVAC 1108 com- 
puter, does not efficiently solve complex systems of equations of the type 
given by equation ( 1 ), (i.e., steady state time harmonic rigid format 8 ) 
when the double precision option is used. Experience on large problems 
(e.g., the size of the one in ref. 2, pg 435) has demonstrated that in the 
solving of the equation (1), NASTRAN has spent 1 iterally. hours in the decom- 
position operation. In order to obtain reasonable run times, a PARAM 
DEC0M0PT4 card is added to the bulk data in order to force the decomposition 
to work in single precision. In comparing the NASTRAN Run 2 (Table 1) 
(1108-S.P. ) results to the exact solution over the frequency range 
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(3. - 8. kHz) where the mesh is sufficiently fine,* it is noted that at 6.0 
kHz, a 154.0% error in the transmitted pressure is experienced. The corre- 
sponding error on the scattered pressure (Run 2, Table 2) is not as severe, 
namely 16.5%. It is noted that the 154% error is the situation that moti- 
vated this entire comparative study. The percent errors at 3. kHz, 4.0 kHz 
and 8 kHz are also much larger than should be expected for the mesh size 
employed. The largest residual, | A | , in the residual vector {A}, is .272 at 
3.0 kHz; this is in comparison to a load vector component of the size (.20) 

(the largest residual was not at a loaded node however). This is a further 
indication that the solution resulting from decomposition in single precision 
is not accurate enough. Upon running the same problem on a Level 17.0 version 
of NASA-NASTRAN on an 1108 computer (single precision must be invoked with 
DMAP) in single precision. Run 3a still results in a similar bad solution with 
a similar worse residual | A | . The fact that Level 17.0 uses a different 
decomposition algorithm did not improve the bad results. However, again run- 
ning the same problem on Level 17.0 on the UNIVAC 1108 in double precision 
(the default situation) resulted in excellent results in comparison to the 
exact solution. For example, at the 6.0 kHz, the percent error reduced from 
154.0% down to 0.4% error. Similarly good results were obtained in the entire 
(3.0 - 8.0 kHz range) for both the transmitted and scattered pressure. The 
double precision gave a very small worse residual at 3,000 Hz, |a| = 2.099xl0“ 8 , 
which suggests that equation (1) has been solved accurately. In comparing CPU 
times between Run 3a and Run 4a in Table 1, it appears there is little penalty 
in CPU time between single and double precision runs for Level 17.0. Conse- 
quently, there is really no incentive (from a time saving point of view) to 
make Sol. 8 type runs in single precision, as there was for Level 15.5 NASTRAN. 

Again running the same problem on Level 48B, MSC version of NASTRAN on a 
CDC CYBER 175 computer (Run 6), very good results were obtained in relation 
to the exact solution over the (3.0 - 8.0 kHz range); further, excellent con- 
sistency with the UNIVAC 1108 double precision runs is demonstrated at all fre- 
quencies by comparing Run 4a with Run 6. The worst residual, |a[, on the CDC 
computer is not as good as the 1108 double precision run, but this is expected 
since the CDC single precision word length is slightly smaller than the 1108 
double precision word length; however, it should be noted that differences in 
the decomposition algorithms could also account for differences in the worse 
residual, even if the word lengths were the same. 

The Level 52, MSC version of NASTRAN on a DEC-VAX computer (Run 7) gave 
comparible results to the Level 17.0 NASA NASTRAN double precision 1108, and to 
Level 48B, MSC NASTRAN on a CDC-CYBER 175. The word length in double precision 
on the VAX is slightly more than the single precision CDC and slightly less 
than the double precision 1108. It is noted the Run 7 gave the smallest worse 
residual . 


* It should be emphasized that in comparing any NASTRAN result to the exact 
solution, for the purpose of measuring the formulation quality, it should be. 
done only in the frequency range of 3.0 - 8.0 kHz where the mesh satisfies the 
10 element/wave length criterion. 
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Computer Type/Level of NASTRAN Sensitivity 


It is inconvenient to run the same exact version of NASTRAN on different 
computers due to leasing restrictions,, consequently this combination was not 
done. Thus running NASTRAN on different computers always involved running a 
different NASTRAN version as well. In scanning the results of Table 1 and 
Table 2, all the runs performed with decomposition precision using word lengths 
between 60 - 72 bits (i.e., Runs 3a, 4a, 6, 7) gave both accurate results in 
comparison to the exact solution (3.0 - 8.0 kHz range) and consistent results 
from machine-to-machine and version-to-version. We have purposely not commented 
on the accuracy of Runs 5, 8 and 9 involving the DEC-VAX computer due to some 
reservations we have regarding the operating conditions of the particular VAX 
on which these runs were made, and will be discussed next. 


Solution Repeatability 


During the process of preparing this collection of comparative runs, an 
unexplained phenomenom (which is still unexplained as of this writing) occurred, 
namely the fact that Level 60 MSC-NASTRAN run on a DEC-VAX computer gave dif- 
ferent results to the same problem upon rerunning the same data. The Run 9 
series of runs were made on rigid format 26 which is comparable (there are dif- 
ferences in the decomposition routine) to NASA-NASTRAN rigid format 8 „ For 
example the input producing Run 9a was resubmitted over again (producing Run 
9b) so that the residual vector {A} could be printed (employing the DMAP 
instructions given earlier). In comparing the solutions, the results were 
slightly different (e.g.,,5.2% at 3.0 kHz). The same input data was again 
rerun on successive days producing Runs 9c and 9d. Run 9c is the closest to 
the more stable results made on the 1108 and CDC computers. 

Using the same DEC-VAX computer facility, the base case input was resub- 
mitted again employing the Level 17.5 NASA/ GODDARD NASTRAN and thus producing 
Runs 5a and 5b. Again the nonrepeatability of the solution on the VAX was 
experienced, this time with an entirely different version of NASTRAN. 

An MSC Level 60 DEC-VAX computer run was made similar to Run 9a through 
Run 9d, except that the older rigid format 8 instead of the newer MSC rigid 
format 26 was employed. The results were poor in comparison to the exact solu- 
tion, further, the worse residual of | A | = .2127 was unacceptably high. No 
reruns of the same input were made on this version and level of NASTRAN. 

In order to demonstrate that repeatable results are possible (a notion we 
usually assume is true on most modern computers), the base case data was rerun 
on the UNIVAC 1108 computer; single precision Runs 3a and 3b were totally 
repeatable as well as double precision Runs 4a and 4b c 

The facts that (1) the VAX computer resulted in nonrepeatable results 
employing two separate versions of NASTRAN and (2) the Level 60 MSC NASTRAN for 
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the VAX computer, Sol. 8 (Run 8) gave poor results, strongly suggests a prob- 
lem with the particular VAX computer on which the runs were made. The follow- 
ing list provides some possible reasons for nonrepeatabil ity: 

• computer central processing drops a bit in the main memory or 
operating register 

• main memory itself drops a bit between storing and retrieving data 

• floating point accelerator drops a bit during calculations 

• disk subsystem (drive or interface drops a bit) 

• computer temperature rises due to air conditioning not keeping 
up with thermal load during the summer months when the runs were 
made; this could result in dropping a bit by one of the four 
above mentioned possibilities. 

Finally it is noted that as of this writing, non NASTRAN users of the same VAX 
computer that made the runs reported here, did not report any repeatability 
problems with computer program results totally unrelated to NASTRAN. 


CONCLUSIONS 


The paper is concerned with both the accuracy of equation solvers and 
with the accuracy of the problem formulation for large dynamic steady state 
problems (e.g., rigid format 8). Based on a series of computer runs on dif- 
ferent versions of NASTRAN on different computers, the following set of con- 
clusions are drawn: 

• NASTRAN solution decomposition algorithms employing less than a 60 
bit word could lead to serious errors in the results (e.g., employ- 
ing 36 bit single precision words on an 1108 computer gave up to 
154% error in the solution with both Level 15.5 and Level 17.0 ver- 
sions of NASA NASTRAN. 

• Correlation between results run on the 1108 double precision Level 
17.0 NASA NASTRAN; C DC-CYBER 175 single precision Level 48B MSC- 
NASTRAN; and DEC-VAX double precision Level 52 MSC-NASTRAN were 
excellent. 

• Unexplained unrepeatability of results were experienced on the DEC- 
VAX computer for both MSC-Level 60 NASTRAN and NASA-G0DDARD Level 
17.5 NASTRAN; Level 17.0 of NASA NASTRAN had no repeatability prob- 
lems on an 1108 computer. 

• A minimum of 10 elements per wave length should be used to model 
traveling wave propagation problems of the type treated in this 
paper; coarser meshes lead to increasingly bad results when comparing 
NASTRAN results to the exact continuum solution to the same problem. 
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As a final comment regarding repeatability, it is noted that it is not being 
suggested that this problem is one to be found in all VAX computers employing 
NASTRAN. The spirit of the NASTRAN Colloquium is to share USER 1 s experiences, 
thus it was felt that our problem should be brought to the attention of the 
NASTRAN USER'S community in the event .that similar problems are encountered by 
others. 
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FIGURE 1 - TEST PROBLEM 



FIGURE 2 - FINITE ELEMENT MODEL 



• 856 Elements (CQDMEM, (TRMEM, CBAR) 

• 588 Grid Points 


element 100378 





SoPo (Single Precision) 
D.Po (Double Precision) 


TABLE 1 - NASTRAN COMPARATIVE SOLUTIONS 
(BACK SIDE PRESSURE, element 100352) 




S.P. (Single Precision) 
D.P. (Double Precision) 


TABLE 2 - NASTRAN COMPARATIVE SOLUTIONS 
(FRONT SIDE PRESSURE, element 100378) 


FREQUENCY 

FRONT SIDE SCATTERED PRESSURE, Ps/Pq (amplitude) 


Dl IM 

SOLUTION SOURCE 

35. kHz 

22.5 kHz 

17.5 kHz 

8. kHz 

6. kHz 

4. kHz 

3. kHz 

NUMBER 

EXACT ANALYTICAL (Ref. 1) 

.04253 

.08580 

.13310 

.34059 

.44266 

.59075 

.59907 

1 

Level 15.5 NASA-NASTRAN, UN I VAC 1008 
Computer, 36 Bit S.P. Decomposition (Sol .8) 

.15093 

.05609 

.11772 

.32372 

.36949 

.59695 

.62678 

2 

Level 17.0 NASA-NASTRAN, UNIVAC 1108 

.15093 

.05603 

.11759 

.32875 

.39233 

.58218 

.60932 

3a 

Computer, 36 Bit S.P. Decomposition (Sol .8) 

.15093 

.05603 

.11759 

.32875 

.39233 

.58218 

.60932 

3b 

Level 17.0 NASA-NASTRAN, UNIVAC 1108 

.15093 

.05592 

.11744 

.33419 

.44176 

.59134 

.59943 

4a 

Computer, 72 Bit D.P. Decomposition (Sol .8) 

.15093 

.05592 

.11744 

.33419 

.44176 

.59134 

.59943 

4b 

Level 17.5 NASA/GODDARD NASTRAN, DEC-VAX 

.14389 

.06016 

.12378 

.29474 

.40715 

.58416 

.59813 

5a 

Computer, 64 Bit D.P. Decomposition (Sol .8) 

.14385 

.06717 

.11805 

.24716 

.32194 

.50502 

.56985 

5b 

Level 48B MSC-NASTRAN CDC CYBER 175 
Computer, 60 Bit S.P. Decomposition (Sol. 8) 

.15093 

.05591 

.11744 

1 

.33419 

.44176 

.59134 

.59942 

6 

Level 52 MSC-NASTRAN, DEC-VAX Computer, 
64 Bit D.P. Decomposition (Sol -8) 

.15093 

.05592 

.11744 

.33419 

.44176 

.59134 

.59942 

7 

Level 60 MSC-NASTRAN, DEC-VAX Computer, 
64 Bit D.P. Decomposition (Sol. 8) 

.16173 

.05591 

.11745 

.33203 

.43973 

.59726 

.59972 

8 

Level 60 MSC-NASTRAN, DEC-VAX Computer, 

.15304 

.05422 

.11353 

.28611 

.44209 

.55800 

.50161 

9a 

64 Bit D.P. Decomposition (Sol. 26) 

.15100 

.05838 

.11587 

.33329 

.44178 

.59291 

.60261 

9b 

(4 runs on same data) 

.15095 

.05592 

.11822 

.33563 

.44131 

,59211 

.59944 

9c 


.15093 

.05592 

.12031 

.33385 

.45400 

.59269 

.59942 

9d 
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ABSTRACT 


The stresses in the CTRAPRG and CTRIARG ring elements are not calculated 
for any of the dynamic solutions in the current COSMIC version of NASTRAN. 

This paper presents a DMAP alter sequence for Solution 8 and post-processing 
program, NASTPOST, to calculate these stresses. Test cases are presented 
which describe the method. The stiffness and the consistent versus concen- 
trated mass problems which have been ascribed to this element are reviewed. 

The DMAP alter sequence introduces Solution 8 displacements to a Solution 
1 module to calculate Real and Imaginary stress components during the execu- 
tion of Solution 8. The post-processor, NASTPOST, calculates the magnitude/ 
phase stress results. 

The DMAP sequence has been written specifically for Level 52 MSC/NASTRAN, 
but can certainly be used for any COSMIC version with slight modification. 


INTRODUCTION 


None of the currently documented versions of NASTRAN calculate the 
dynamic stresses in the CTRAPRG and CTRIARG solid of revolution elements. The 
stresses for these elements are calculated in NASTRAN for static solutions 
(e.g.. Solution 1) but not in the dynamic solutions (e.g., Solution 8)„ Com- 
ments have been made by others which express the reasons for not including the 
stress calculations are related to the formulation of the mass matrix for the 
element. 

Sample problems are given to show that the difference between the consist- 
ent and concentrated mass approach is greater than one might expect from argu- 
ments solely between the merits of consistent or concentrated mass. 

This paper describes a DMAP alter sequence for Solution 8 and a post- 
processing program, NASTPOST, to calculate these dynamic stresses. The DMAP 
alter sequence introduces the displacements computed in Solution 8 to a Solu- 
tion 1 module to calculate the complex stresses in the form of real and 
imaginary components. The post-processor, NASTPOST, calculates the stresses 
in the form of magnitude/phase. 
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DISCUSSION 


It is not spelled out in the NASTRAN Users Manual that stresses for the 
solid of revolution elements are not calculated for dynamic solutions. There- 
fore, if one asks for stresses in a Solution 8 case control, the run is not 
aborted, but no stresses are obtained. 

In order to perform noise path studies of an axisymmetric structure it 
became necessary to obtain these stresses. At first, the displacements for 
the entire structure, obtained from a Solution 8 forced vibration analysis 
were written into an output file; then these displacements, less one, were 
written into SPC format as enforced displacements for a Solution static analy- 
sis (this was done for the real and imaginary components separately). This 
technique was later modified, utilizing the DMAP alter sequence A0S8$CS and a 
post-processor, NASTPOST. 

The DMAP alter sequence is given in Figure 1. The major points are: 

• The user can specify output requests as usual for SPCFORCES 
and DISPLACEMENTS. 

• The user should specify STRESS (PUNCH) = ALL or a particular 
set ID if he wishes to subsequently use NASTPOST to calculate 
the magnitude/phase. This punched file will be sent to the 
users system space. (FOR 013. DAT for the MSC/NASTRAN VAX 11/780 
VERSION). 

• A0S8$CS should be pi deed on the user's RFALTER library and 
executed then by calling RFAI = A0S8$CS. 

The program NASTPOST is given in the appendix and is used to calculate 
magnitude/phase stress components from real/imaginary stress components. The 
major points are: 

• The components from F0R013.DAT above, are used as input to 
calculate the magnitude/phase stress components. 

• This program can be run immediately after the execution of 
MSC/NASTRAN or at some later time. 

The test problem for A0S8$CS and NASTPOST is a circular plate fixed at 
the edges and driven by a single force, 100 dynes, at the center, normal to 
the plane of the plate. The finite element control model is the CQUAD2 and 
CTRIAG2 bending element model shown in Figure 2. The CTRAPRG model, shown in 
Figure 3, is formulated as a concentrated or consistent mass for each of the 
runs. The NASTRAN default value is the consistent mass matrix. The concen- 
trated mass matrix is entered as C0NM2 data. The three cases are compared in 
Table 1 for static, 2000 Hz and 8000 Hz at a position near the concentrated 
load and at the fixed edge. 
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The concentrated mass formulation gives good results, as compared to the 
control model. The consistent mass, or default formulation, gives results 
which do not agree with the control model at either the low, 2 kHz, or high, 

8 kHz, forcing frequencies. 

The static solution agrees very well with the control model which indi- 
cates that the stiffness of the model is represented correctly by solid of 
revolution elements. The error therefore is associated with the mass matrix 
formulation. The degree of error is obviously greater than one would expect 
from the normal arguments of consistent versus concentrated mass differences. 1 

It can be argued that the use of cyclic symmetry with 3D elements rather 
than solid of revolution elements would have been a possible solution. This 
is certainly an avenue that deserves added investigation for comparison of 
cost and accuracy of solution compared to the solid of revolution elements with 
concentrated mass matrix. 


CONCLUDING REMARKS 


A DMAP alter sequence for Solution 8 and a post-processing program 
NASTPOST has been presented to calculate the dynamic stresses in CTRAPRG and 
CTRIARG solid of revolution ring finite elements. Users of this technique are 
cautioned to use the concentrated or lumped mass matrix rather than the con- 
sistent mass (default value) matrix. 

The DMAP sequence has been written specifically for Level 52 MSC/NASTRAN, 
but can certainly be used for any COSMIC version with slight modification. 
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TABLE 1 


COMPARISON OF STRESSES, 3/8 cm from CONCENTRATED LOAD 


FREQUENCY 

01 

2 kHz 

8 

QUAD2 

134.4 

75.5 

66.4 

TRAPRG (CONS.) 

132.3 

17.2 

63.1 

TRARG (CONC.) 

132.3 

96. 

60.5 


TABLE 2 

COMPARISON OF STRESSES, 3/8 cm from FIXED EDGE 


FREQUENCY 

O 1 

2 kHz 

8 

QUAD2 

44.4 

34.2 

38.2 

TRAPRG (CONS.) 

45.6 

27.0 

10.0 

TRAPRG (CONC.) 

45.6 

33.0 

36.0 


1 OBTAINED FROM SOLUTION 1 
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FIGURE 1 - ALTER A0S8$CS 


BEGINNING OF ALTER A0S8ICS 

THIS ALTER PACKAGE IS USED TO CALCULATE 

^DISPLACEMENTS (REAL/ 1 MAG I NARY ) OR 

(MAGNITUDE/PHASE) 

*SPCFORCES (REAL/IMAGINARY) OR 

(MAGNITUDE/PHASE) 

^STRESSES (REAL/IMAGINARY) 

FOR THE CTRAPRG AND CTRIARG RING ELEMENTS 


CASE CONTROL INPUT 



FIGURE 1 - (Cont'd) 


THE USER SHOULD SELECT THE DESIRED 
OUTPUT AS USUAL FOR DISPLACEMENTS 
AND SPCFORCES. 

THE USER SHOULD SELECT THE PUNCH 
OPTION FOR STRESS IF IT IS DESIRED TO 
SUBSEQUENTLY CALCULATE (MAGNITUDE/ 
PHASE) USING A POST-PROCESSING PROGRAM 


ALTER 166 

OFP OPPC1 , OQPC 1 , OUPUC1 , , , //U # N # CARDNO « 

ALTER 185,186 
PARAM //STSR/13/-64 f 

GP3 GE0M3, EQEXIN, GE0M2/ , ETT/O/U, N, NOGRAU/0 % 



FIGURE 1 - (Cont'd) 


cn 

i£> 


PARAHL UP0C//C,N,TRAILER/2/M»N,R0US * 

PIATGEN ,/UNIT/l/ROUS * 


MODTRL UPUC////3 * 

HPYAD UNIT,UPOC,/ASQR/ * 
DIAGONAL ASQR/ATRN// * 

ADD UPUC,/BSQRA0.8»-1*O> * 


DIAGONAL BSOR/BTRH// * 

SDRS CASECC,CSTH,NPT,DIT,EQEXIN,SIL,ETT,EDT,BGPDT,,, 

XYCDB/, » ,0ESCR, , /STATICS/S, N,N0S0RT2 * 

SDR2 CASECC, CSTH, HPT , DIT, EQEXIN , SIL, ETT , EDT , BGPDT , , , 
XYCDB/, , .OESCI „ /STAT ICS/S , N , N0S0RT2 * 

OFP ,.,OESCR,,//S,N,CARDNO t 
OFP ,,, OESCI, .//S.N.CARDNO S 


ATRfl,EST 

BTRH.EST 


» 


» 


PARAH //STSR/7/-64 f 


ENDALTER % 

t 





XXX 



APPENDIX A 


THE NASTPOST PROGRAM 
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c DATA SET NASTPOST AT LEVEL 11/05 79 

COIWON /HDRCOH/TITLEt lfi l.SU8T( 16)»LABEL( 16^ 

DATA DT I T/ ' *T IT ' / , C ASE/ ' C ASE V . DS J»/ *SUt / . 
t DELE/ ' 0ELE ' /. DSTR/ ' STR'/.DLAD/ *LA» / 

DATA 1036. 1037/2**/ 

1 CONTINUE 
REUIND 7 

C - GET TITLE CARD 

5 CONTINUE r 

RSA P<7, 960,ENP»999> TEMP. TITLE 
IF(TENP.EO.DTIT) GO TO 6 
GO TO 5 

C - GET SUBTITLE CARD 

* READ<7^9— ,END*999) TENP.SU1T 
IFtTENP.EO.SSUS) GO TO 7 
GO TO 6 

C - GET LABEL CARD 

7 re5Si7%.ehd-m») temp.lasel 

IF(TEA*.EQ.DLA8) GO TO 1* 

GO TO 7 

C - GET STRESS CARS 

1* CONTINUE ___ 

REAP(7.91*. END-9991 TENS 
IFCTENP.EG.DSTR) GO TO 2* 

GO TO 10 

C - GET SUSCASE IDENTIFICATION 

REAsIfflsO.ENP-O**) TEMP.ISIS 
IFdEMP.EG.CASE) GO TO 3* 

GO TO 2* 

C - GET ELEMENT TYPE 

IF ( TEMP. NE. DELE) GO TO S 
C - CHECK ELEMENT TYPES 

IFdELTYP.EG.3S) GO TO 3S* 

IF(IELTYP.EG.37» GO TO 37* 

00 TO S _ 

C - ELEMENT TYPE • 3S 

** trttflflgF.M- *) CALL RU3SdSIS.IELTVP.IE0F) 

IF(I03S .EG. I) CALL RC36dSID.IELTVP.IE 0F ) 

IF (1036 .EG. 1 .AMS. I EOF .EG. 1) GO TO 99* 

I03S - H0SCIO3S+I.2) 

GO TO • _ 

C - ELEMENT TYPE • 37 

379 TTtlttrf .EG. A) CALL RU37dSID.IELTYP.IE0F) 

1FCI037 .EG. I) CALL RC37( ISIS. IELTyP t IEOF ) 

IF1I037 .EG. 1 .AMO. IEOF .EG. 1) GO TO 90S 
1037 - N00d037+1»2) 

GO TO S 


000 FOAMATt A4.SX.1SA4. At) 
91* F0ANAT(8X.A4) 

980 FORMAT! 4X.A4.8X, 19) 
939 


10001 

00002 

0O003 

0000-4 

00005 

00006 
000*7 


00009 

0001 * 

00011 

00*12 

00*13 

00014 

0*015 

00016 

00*17 

0*018 

00019 

00*2* 

00021 

00022 

00*23 

00024 

00*25 

00026 

00027 

00*28 

00*29 

00030 

0*031 

0*032 

00033 

00*34 

00035 

00036 
0*037 
00*38 
00*39 


0*041 

00*42 

00*43 

00*44 

0*045 


00*47 

0*048 

0*049 


00051 

0*052 

0*063 


***5S 

00056 

00057 
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OOOS 9 


FORMAT! A4, 12X, 111 ) 

END 

C DATA SET NASTRU36 AT LEUEL 004 AS OF 11/02/79 

SUBROUTINE RU36! ISID, IELTVP, IEOF ) 

DIMENSION TEMP!2)»DATA!4) 

DATA TITLE/'BT ' / , CONT/ ' -CON ' / , BLANK/ ' '/ 

DATA INN, I OUT/?, 9/ 

REWIND I OUT 
PRINT 10 

10 FORMAT! 'SUBROUTINE RU36' ) 

READ! INN, 900, END-999 ) IELNO.DATA! 1 ),DATA!2 >. DATA! 3) 

#01 CONTINUE 

READ! INN, 910, END-990 ) CARDN, DATA! 4) 

IF (CARDN .NE. CONT) CO TO 990 
URITE! IOUT ) ISID, IELTYP, IELNO.DATA 
C READ! INN, 920, END-999) TEMP 

C BACKSPACE INN 

CALL BACKSP ( TENP, INN, 1999 ) 

IF !T£HP(l).EO. BLANK) 

S READ! 10, 900. END-999 )IELNO, DATA! 1 ), DATA! 2), DATA! 3 ) 
IF!TENP!1) .EO. BLANK) 60 TO 001 
IF1TENP11) -NE. TITLE) CO TO 990 
800 CONTINUE 

ENDFILE IOUT 
REMIND IOUT 
RETURN 

990 CONTINUE 
STOP 3600 
999 IEOF • 1 
GO TO 800 

900 FORNAT!I10.8X,3£18.6) 

910 F0RNAT(A4.14X.3£18.S) 

920 F0RNAT(2A2) 

END 

C DATA SET NASTRU37 AT LEUEL 004 AS OF 11/02/79 

SUBROUTINE RU37! ISID, IELTVP. IEOF) 

10 FORMAT! 'SUBROUTINE RU37' ) 

DIMENSION TENP!2).DATA!20),KKREAD!33> 

DATA TITLE/ '9T '/.CONT/' -CON'/, BLANK/' '/ 

DATA INN, IOUT/7,8/ 

REWIND IOUT 
PRINT 10 

READ! INN, 900. END-999 ) IELNO, DATA ! 1 >, DATA ! 2 ) , DATA! 3 ) 

001 CONTINUE 

READ! INN, 910, END-990) CARDN, DATA!4 ), DATA! S),DATA!6) 

IF! CARDN .NE. CONT) CO TO 990 

READ!INN,910,END>990) CARDN. DATA! 7), DATA! 8). DATA! 9) 

IF! CARDN .NE. CONT) GO TO 990 

READ! INN. 910, END-990 ) CARDN, DATA ! 10 > , DATA! 11), DATA ! 12 ) 
IF (CARDN .HE. CONT) 00 TO 990 

READ! INN, 910.END-990) CARDN, DATA! 13), DATA! 14 ), DATA! IS ) 
IF (CARDN .NE. CONT) GO TO 990 

READ! INN, 910, END-990 ) CARDN, DATA! IS ), DATA! 17 >, DATA! 18 ) 

IF (CARDN .NE. CONT) GO TO 990 

READ! INN, 910, END-990) CARDN, DATA! 19 ),DATA(20) 

IF (CARDN .NE. CONT) GO TO 990 
UNITE! IOUT) IBID, IELTVP, IELNO, DATA 
C READ! INN, 920, END-999) TENP 

C BACKSPACE INN 

READ! INN, 930,END-999)KKREAD 
REWIND 10 

WRITE! 10,930 IKKREAD 
REWIND 


• -S 


00001 

00002 
00003 
30004 
00005 


00007 

00008 

00009 

00010 
00011 
00012 


00013 

00014 

00015 

00016 

00017 

00018 

00019 

00020 
00021 
00022 

00023 

00024 

00025 

00026 

00001 


00003 

00004 

00005 

00007 

00008 

00009 

00010 
00011 
00012 

00013 

00014 

00015 

00016 

00017 

00018 

00019 

00020 
00021 
00022 
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1* 


R£AD(10.92O)T£nP 
REWIND 1* BI 

» X READUo! 9 w!eAd® 999^ I E LMO. DAT A( 1 ) . DATA ( 2 > , DATA ( 3 > 
IF(TEHP(1>.EQ. BLANK) 6GT0 Ml 
IF<TEMP(1> .HE. TITLE) ^0 TO 990 
800 CONTINUE 

ENDFILE I OUT 
REWIND I OUT 
RETURN 
99* CONTINUE 
STOP 3700 
999 I EOF • 1 
GO TO 800 

g«0 FORIIATI 1 1#»8X» 3E18.fi) 

910 F0RNAT<A4,14X,3E18.6> 

980 F0RNATC2A2) 

930 FORMAT (33A4 ) 

• SUIROUT INE^RCsSdS ID?IELT VP ,VloF ) « “ " 

“ 01 ‘ cl^ C 4^ RHAG 14). PHASE C > 

am isasKio^*?:^ ' 

PRINT 10 
IELCNT * 99 

KAD< INN.U#!InD?99«) IELNO.DATAI < 1 > . DATAI ( 8 ) . DATAI ( 3 ) 

•“ R€ADnNN?910.END-990 ) CARDN^DATAI (4 ) 

«O T0%9# 

IF(IELNOR .NE. IELNO) CO TO 990 
° 0 RAA6<i> - 1 SORT< DAT AR (I ’ 5 

K'MKii! -a- S:S! a. T sl,!r 


DATAI(I>*DATAI(I)> 


IF (DATAI (I ) .EG. W.w> Lnss;;; 

IF(DATAKI) .CT. 0.0 PHASE 
IF (DATAI (I) .LT. #.§) PHASE(I) 

GO TO «98 

Satio^abscdataki VDATARCI)) 

X IF(PATAI< I )‘l T. 0.0 .AND. ,PPTAR ( 1 l.LT.O.O) 

X lF(DWAICI)’l.2S IT*' ®JTAR«).8T.#.G) 

X PHASC(l) • PHAfiE(I) ♦ 870.0 

c 699 SmoPRT.93G) wi&ielweuio.datar.datai 

IF (IELCNT .LT. SO) 00 TO TOO 
CALL HKMOSIO) 

IELCNT • 0 

"• tssn «l»t . . 


0.0 

90.0 

870.0 


<*0054 

00025 

00026 

30027 

00028 

00029 

00030 

00031 

00032 

00033 

00034 

00035 

00036 
00001 

00008 

00003 

00004 

00005 

00006 
00008 

00009 

00010 
00011 
00012 

00013 

00014 

00015 

00016 

00017 

00018 

00019 

00020 
00021 
00028 
00023 
00084 

00025 

00026 

00027 

00028 


00031 

00038 

00033 

00034 

00035 

00036 
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94437 

URITEUPRT.940) «WO.<^JJ« <I, « p H ASC,,,, ' I ‘ 1 ' 4> 
READ(INN,929,END*999> TENP 

“IS'ilcJsPCTEW. 

IF(TEHPtl) .EO. BLANK) CO TO Ml 
IF(TEf1P< 1 ) .HE* TITLE) GO TO 99* 

RETURN 
999 CONTINUE 
STOP 3691 
999 IEOF • 1 

dm FORtMT<I19*8X,3£18.6) 

91# F0RmT(A4.14X.3£18.6> 

> Zmvssa&x* ™ “ ' li '* 5 ' 79 

MTA IPRT,INN,I0UT/«.7.8/ 

PRINT 19 
IELCNT ■ 19 

REAof INN > 999?END^999 ) IELNO, DATAI ( 1 > * DATAI <2 ) > DATAI ( 3 ) 

• #1 re^^m.en^ > c^dn^ata:< 4 > * DATftI 15 > • DftTftl (6 ‘ 

RE ADdNN, 919 ", END^WO^C ARDN?DATA I* 7 ) , D AT A I ( 8 ) # D AT A I ( 9 > 

READHNN. 919:END^999^CAR0N?SaTAI< 13 ) , DATAI ( 14 ), DATAI 
RErS?INN, 919:ENB^999^ C ARD^DATAI 1 18 ) # DATAI C 17 ) » 9ATAI C 1 8 ) 
^^,9f.:^K^TAm9>.»ATAn29. 
5^S5 t-, , 1siMlMe^r.i»atar 

IF(ISIB *NE. ISIDRIGO TO 999 
IFCIELTYP .HE. IELTPR1 CO T0999 
1FUELN0R .HE. IELNO) CO TO 999 

R» :§: W £${£ u\ 

S!mt*!<!> :»: W! WBil! • 

00 TO 889 

RATIO^AtSIDATAI < I )/DATAR«I ) ) 


DATAICI )*DATAI<I>) 


03038 

0O039 

00040 


03041 

33042 

00043 

00044 
300 45 
30046 

00047 

00048 

00049 

00050 

00051 

00052 

00053 

00001 

00002 

00003 

00004 

00005 

00006 
00008 

00009 

00010 
00011 
00012 

00013 

00014 

00015 

00016 

00017 

00018 
00019 
00029 
00021 
00022 

00023 

00024 

00025 

00026 
00027 

00928 

00929 
00939 

00031 

00032 
00933 
0M34 

00035 

00036 

00037 
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X PHASE ( I ) • PHASE! I ) + 180.0 

IF t DATA I ( I ) .LT.9.9 .AMD. DATARtl ).GT.0.0) 

X PHASE ( I ) • PHASE ! I ) + 27*. • 

699 CONTINUE 

: URITE! IPRT, 93® ) ISID.IELTYP, IELNO, DATAR, DATAI 

IF< IELCNT .LE. 7) GO TO 790 
CALL HD37CISID) 

IELCNT • 9 
799 CONTINUE 

IELCNT • IELCNT ♦ l 
DO 710 I • l.S 
J • 4*!I-1> ♦ 1 
K • J + 3 

IF< I .EO. 1> URITE! IPRT, 940) IELNO.I. 

X ( <RNAG( 1X1 ), PHASE! IX 1 ) )>IX1 ■J.K) 

IF( I .NE. 1) URITE!IPRT,959) I, 

X ! (RMAG! 1X1 ), PHASE! 1X1 )),IX1*J,K) 

719 CONTINUE 

WRITE! IPRT, 960) 

5 READ! INN. 929. END-999) TEMP 

: BACKSPACE INN 

CALL BACKSP! TEMP, INN, 1999 ) 

9 READ^lS! 990^END*999 ) I ELMO, DATAI 11), DATA I ( 2 ) , DATAI C 3 ) 
IF1TENP! 1 ) .EO. BLANK) 00 TO 991 
IF(TENP<1 ) .NE. TITLE) 00 TO 999 
RETURN 

999 CONTINUE 
STOP 3791 
999 IEOF « 1 
RETURN 

999 F0RHAT!I19,8X.3E18.6> 

919 F0RNAT(A4, 14X.3E18.6 ) 
q m F0fiMT(2A2) 

939 FORMAT! IX, 31 10, 10!/. 4! 5X, 1PE13.6) ) ) 

949 F0RNAT!1X.IS.1X,I3.4X,4!1PE12.5,' /' .8PF19.S.5X ) ) 

9S9 F0RMAT(7X.I3.4X.4!1PE12.S.' /',0PF10.S.5X> > 

969 FORMAT! ' ' ) 

END 

C DATA SET NASTHD36 AT LEUEL 997 AS OF 19/24/79 

SUBROUTINE HD361ISID) 

19 FORMAT! 'SUBROUTINE HD36' ) __ _ 

COMMON /HDRC0M/TITLE!IS)»SUBT(16),LABEL(16> 

PRINT 19 
IPRT • 6 

WRITE 1IPRT. 199) TITLE 
WRITE! IPRT, 119) SUBT 
WRITE! IPRT, 129) LABEL, IBID 
WRITE! IPRT. 149) 

WRITE! IPRT, 150) 

WRITE! IPRT. 169) 

WRITE! IPRT, 179) 

RETURN 

199 FORMAT! '1',3X,1SA4,AB> 

119 FORMAT!' ',3X,1SA4,A2) _ 

129 FORMAT! '0',3X, 15A4,A2,S0X, 'SUBCASE'. 13) 

FORMT( * * ) 

149 F0RMAT(27X, 'STRESSES FOR THE TRIAN 


03338 

03339 
aoo«0 

30041 

00042 

00343 

30044 

0334S 

30046 

00047 

00048 

00049 

00050 

00051 
30052 

00053 

00054 

00055 
00956 

00057 

00058 


00059 

00060 
30061 
00062 

00063 

00064 

00065 

00066 

00067 

00068 
00069 
00979 
09971 

00972 

00973 

00091 

00002 

90993 

00994 

00995 
90996 

00997 

00998 
00099 
00019 

00911 

00912 
00013 

00914 

00915 

00916 
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X 'G U L A R RINGS ICTRIAR6)') 

159 F0RMATI61X, ' !NAGNITUDE/PHAS£)' ) 

169 FORMAT ( -IX. ' EL ' , 17X. 'RADIAL ' . 19X, 'CIRCUMFERENTIAL ' , 

X 19X, 'AXIAL', 24X, 'SHEAR' > 

179 FORMAT! 4X, 'ID' . 19X. ' !X >' ,24X, '! THETA )',24X, ' (2 ) '. 

X 26X, ' (ZX)' ) 

END 

C DATA SET NASTHD37 AT LEVEL 996 AS OF 19/24/79 

SUBROUTINE HD37USID) 

19 FORMAT! 'SUDROUTINE HD37' ) 

COMMON /HDRCOM/TITLE 116), SUIT! 16), LABEL! 16 ) 

PRINT 19 
1PRT • 6 

URITE! IPRT. 199 ) TITLE 
URITEC IPRT, 1 19 ) SUIT 
URITE! IPRT. 129) LASEL.IS1D 
WRITE! IPRT. 149) 

WRITE! IPRT. 1S9) 

WRITE! IPRT. 169) 

URITE! IPRT. 179) 

RETURN 

199 FORMAT! '1',3X.1SA4. A2) 

119 FORMAT!' '.3X.1SA4.A2) 

129 FORMAT! '9'.3X. 1SA4.A2.59X. 'SUBCASE'. 13) 

13# FOAfMT ( ' # ) 

149 F0RMAT!27X,'S TRESSES FOR THE TRAPE', 
X 'ZOIDAL RINGS ICTRAPRG)') 

1S9 F0RHAT!61X, '! MAGNITUDE/PHASE ) ' ) 

169 FORMAT! 4X. 'EL' ,2X. 'ST' .13X, 'RADIAL' , 19X, 'CIRCUMFERENTIAL' . 
X 19X. ' AXIAL'. 24X, 'SHEAR') 

179 FORMAT! 4X. ' ID' . 2X, 'PT' , 15X, ' ! X ) ' . 24X, ' ( THETA ) ' , 24X, ' !Z ) ' , 
X 26X.MZX)') 

END 

SUBROUTINE BACKSP! TEMP. INN,* ) 

DIMENSION KKREAD!33),TEHP!2) 

READ! INN,930,END a 999)KKREAD 
REWIND 19 

URITE! 19,939 )KKREAD 
REWIND 19 
READ! 19.929 >TENP 
REWIND 19 

939 FORMAT (33A4) 

989 FORMAT !2A2) 

RETURN 

999 RETURN 1 

ENR 


09917 

99*19 

99919 

99*29 

09021 

90922 

09923 

00091 

00992 

00993 
00004 
0009S 

00996 

00997 
00098 
00999 
09910 
00011 
00012 
09013 
00914 

00015 

00016 

00017 

00018 
00019 
00029 
00921 
00022 
99923 



AN ENHANCEMENT OF NASTRAN FOR THE SEISMIC ANALYSIS OF STRUCTURES 

John W. Burroughs 

Civil Design Department, Ontario Hydro 


SUMMARY 


New modules, bulk data cards and DMAP sequence have been added to 
NASTRAN to aid in the seismic analysis of structures. These allow input 
consisting of acceleration time histories and result in the generation of 
acceleration floor response spectra. The resulting system contains, numerous 
user convenience features, as well as being reasonably efficient. 


INTRODUCTION 


At ONTARIO HYDRO, the. primary analysis tool is NASTRAN. This use began 
with the purchase of level 15.1 and continued with 15-5 and SPERRY/NASTRAN. 
Currently MSC/NASTRAN is being implemented as levels 16 and above are not 
available in CANADA. To perform a seismic analysis of nuclear power plant 
structures, the NASTRAN normal modes analysis has been utilized in 
conjunction with two post processors written at ONTARIO HYDRO. The one 
performs a time history method analysis to generate the desired floor 
response spectra, the other performs a response spectrum method analysis by 
utilizing these floor response spectra. While the response spectrum 
processor has stood the test of time, the time history method processor has 
not . 


The time history analysis post processor- was initially conceived to be 
used with simple stick type lumped mass models. However, as the complexity 
of the analyses increased, this processor was unable to satisfy all 
requirements. In addition, the cost of analysis for the more complex 
structures became excessive. About this time, the need to treat problems 
which were subjected to multiple-support excitations became a requirement. 

Considering the shortcomings of the existing post processor, as well as 
future requirements, it was decided to develop an entirely new capability and 
to include it within NASTRAN. This project was then divided into two main 
development stages. The first, which is described in this paper, is a 
replacement for the existing post processor and provides the ability to 
handle any size problem, efficiency, improved output and user convenience 
features. The second stage, development of which is underway, extends the 
first stage to allow for the consideration of multiple-support excitation • 
problems. 
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SYMBOLS 


a(t) 

b j 

Fi(t) 

Se 

k J 

*1 

m j 

m i 

Pi(t) 

0 

s 


acceleration time history at the ground (g). 

structural damping ratio associated with a particular mode 

j- 

generated force time history acting at a specific degree 
of freedom. 

structural damping ratio associated with the 1th element. 

stiffness matrix associated with a specific mode j. 

stiffness matrix associated with a particular element 1. 

modal mass matrix associated with the jtb mode. 

total structural mass associated with grid point i. 

absolute acceleration time history used as a load for the 
response spectra generation. 

equipment damping ratio for which spectra are required, 
displacement components in modal co-ordinates, 
mode shape associated with the jth mode. 

input and modal frequencies at which the response spectra 
will be computed. 


THEORY 


Since the input for seismic activity is usually available as 
acceleration time histories, and force time histories are required, a 
conversion must be performed. Several techniques are available, the one 
chosen here replaces the acceleration time history by a set of force time 
histories according to the equation 


Fi(t) = a(t) (1) 

This results in a force time history at each free degree of freedom 
corresponding in direction to the input acceleration. 

Once the force time histories have been created, a modal transient 
analysis is performed. The resulting relative acceleration time histories 


sa 



are then converted to absolute accelerations prior to the generation of the 
floor response spectra. This is done by adding the input acceleration to 
each of the computed acceleration time histories in the corresponding 
direction. 

Floor response spectra are computed by performing a transient analysis 
for each of a set of one degree of freedom oscillators. The transient 
analysis performed utilizes the solution of separate second order 
differential equations of the following form. 

£i + 2/3co 0 + u> 0 z = mi ) (2) 


<*>0 0 = 


bi 

2mi 


(3) 


w 


i _ 


Ki 

mi 


(4) 


The frequencies o> 0 utilized are a combination of user input frequencies 
and the frequencies determined for the structure. As results are required 
only in modal coordinates, the value of m^ is arbitrarily set- to 1.0. This 
requires then only the solution of the equation 


i i + 2/3o> 0 £ i + w 0 * £ i = P i (t ) ( 5 ) 

The equation of motion, therefore, corresponds to that of a single degree of 
freedom system having the prescribed damping and frequency properties and 
subjected to the prescribed degree of freedom acceleration. The relative 
acceleration is obtained at each time step as follows: 

£i r n+ 1 = 1 -2/3w 0 £; tn + | _ w i n + I (6) 

To obtain the desired absolute accelerations, the computed acceleration is 
added to the above relative acceleration. The maximum acceleration is then 
retained over all time steps. This procedure is repeated for each of the 
designated frequencies and equipment damping input*. The resulting table of 
maximum acceleration versus frequency is the desired floor response spectra. 


Damping must be included in the analysis. Here, the user may specify 
modal damping, uniform structural damping or element structural damping. The 
preferred technique is to use element structural damping. In this case, the 
composite modal damping values will be computed. These values are based upon 
the fraction of the strain energy sustained by each element in the model. 

The modal damping for the jth mode is computed as follows: 
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(7) 


bj = Kt \ty]\ gl 

i<T 

APPROACH 


This facility has been implemented in NASTRAN by means of a DMAP 
program. This program is a modification of that used for Modal Transient 
Analysis (Rigid Format 12). The general problem flow is as shown in 
Figure 1. 

To implement this facility, four new modules were written. Two of which 
precede the transient analysis module (TRD) and two follow. In addition, 
extensive use is made of existing NASTRAN modules in the solution. 

The standard NASTRAN approach is followed for the matrix generation and 
eigenvalue extraction phase. Following this, the equivalent force time 
histories are created. The input acceleration time histories may come either 
from TABLED1 cards, or from a user file where the tables have been 
prestored. A modal transient analysis is then performed. The output from 
this stage consists of relative acceleration time histories. An added module 
will convert these into absolute accelerations. This matrix is subsequently 
transposed and from it the floor response spectra are generated. Finally, XY 
plots of the spectra may be produced. All normal NASTRAN output is 
available, in addition to the output produced by these new modules. 


NASTRAN IMPLEMENTATION 


To implement this capability in NASTRAN, two bulk data cards, four 
modules and a complete DMAP sequence have been developed. The new bulk data 
cards are SDATA and SET1. 

The SDATA card, used to define the input loading and optionally to 
select the data required to generate the floor response spectra, is 
illustrated in the appendix. The SDATA card is selected by the DLOAD case 
control card. If the acceleration is to be combined with other acceleration 
or force time histories, then the DLOAD bulk data card may be used to combine 
them. Each SDATA card may select acceleration time histories for up to six 
degrees of freedom at any one grid point. In addition, data may be provided 
for the generation of floor response spectra. This data includes the 
equipment damping set and the set of grid points at which spectra are 
desired. Miscellaneous data for the control of the analysis may also be 
provided . 

Hie SET! card is used to define the grid points at which spectra are 
computed. It is selected by the SDATA card and is required only if floor 
response spectra are to be generated. The card format is shown in the 
appendix. 
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The acceleration time histories required may be supplied either as 
TABLED1 cards or from prestored tables on a user file. The latter technique 
is preferred when a standard set of time histories is available. 

Four modules were created for this enhancement. They are SCNTL, SAPF, 


STHGNMX, 

and 

SFRG. 

SCNTL 

— 

This module verifies all data input on SDATA cards and ensures 
that the required sets and data tables are available before 
proceeding with the analysis. 

SAPF 

— 

This module accepts the input acceleration time history and 
generates the required force time histories. This is done by 
creating new forms of the DLT and DIT tables. 

STHGNMX 


This module accepts the relative time histories output from TRD 
and creates the absolute acceleration time histories. In 
addition, this module can, at user request, reduce the size of 
the output matrices prior to subsequent output requests. 

SFRG 


This module accepts the transposed absolute acceleration matrix 
and creates the required floor response spectra. This includes 
the generation of data for XY plots as well as printed output. 
This information is then passed to the XY PLOT and OFP 
modules. 

In addition to these bulk data cards and modules, the DMAP sequence 
contains a number of parameters which may be used by the engineer to select 
optional processing paths. In general, the engineer need not use any of 
these parameters as the default values will select the most appropriate 
options. 

All the 
This includes 
and all which 
be printed or 

normal output from a Modal Transient analysis may be requested, 
both relative and absolute accelerations. The primary output, 
is normally required, is the floor response spectra which may 
plotted as desired. 


SAMPLE PROBLEM 


To illustrate the ease with which this enhancement may be used, the 
structure shown in Figure 3 is analyzed and floor response spectra generated 
at grid point 6. As a point of comparison, the floor response spectra 
generated by the previous' post processor is shown in Figure 4. That portion 
of the NASTRAN input data required for this enhancement is shown in 
Figures 5, 6 and 7 for the Executive control, Case control and Bulk data 
decks respectively. The plot of the resulting floor response spectra is 
shown in Figure 8 and a portion of the printed output in Figure 9* A review 
of this data and the resulting output demonstrates the ease with which floor 
response spectra may be generated. 
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EXTENSIONS 


Several extensions to. this enhancement are either planned or in 
progress. One will be to implement this facility in the other NASTRAN 
versions available at ONTARIO HYDRO. Another is the extension to solve 
multiple-support excitation problems. Upon completion of this last item, the 
original reason for. developing these features will be satisfied. All future 
work will then add more user convenience features or efficiency 
improvements . 


CONCLUSIONS ' 


As a result of this enhancement of NASTRAN, an easy to use capability 
for the generation of floor response spectra has been made available. This 
is an extension of the existing Modal Transient analysis which retains all 
the original capabilities. In addition, the approach used is relatively 
efficient in terms of computer resources and user interaction required. It 
is a vast improvement over the original post processor and removes all of the 
restrictions inherent in it. 
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APPENDIX 


Input data card SDATA 

Description: Defines an input acceleration time dependent loading and 

various parameters for the generation of floor response 
spectra. 

Format and example 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

SDATA 

SID 

RUN 

XQT 

DTM 

SAVE 

SOURCE 

THN1 

DIR1 

+abc 

SDATA 

101 

-1 

3 

2 



BRTEQH 

2 

+SD1 


+abc 

IDFSP 

IDEQ 

NSKO 

IDUMP 

THN2 

DIR2 

THN3 

DIR3 

+def 

P| 

10 

20 


. . . . 







+def 


THN4 

DIR4 

THN5 

DIR5 

THN6 

DIR6 

H 













Field 

Contents 


SID 

Set identification number (integer 

> 0) 

RUN 

Run type control parameter < 0 - 

> 0 - 

simple structure 
multi excitation 

XQT 

Erection phase control 



= 1 generate time histories only 
= 2 generate floor response spectra only 
= 3 both 1 and 2 


DTM Identification number of the grid point at which the 

acceleration time history is applied (integer > 0) 
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SAVE 


A flag indicating that the output time histories are to be 
saved for subsequent use (integer > 0 or blank) 

SOURCE Identifier of the file on which the input time history is 

stored (integer ^ 0, default = 0) 

THNi Name of the acceleration time history if on a file, or the 

id of a TABLED1 card. 


DIRi Direction associated with this' time history 

(1 ^ integer 6) 

IDFSP Set id of a SET1 card on which the points at which floor 

spectra as desired are listed (integer > 0 or blank) 


IDEQ 

NSKO 

IDUMP 


Remarks : 


Set id of FREQ) FREQ1 or FREQ2 cards on which the 
equipment dampings are defined (integer > 0 or blank) 

Alternate skip factor to reduce the output time histories 
by (integer > 0 or blank) 

Intermediate output flag (integer > 0 or blank) 

= 1 print input tables 
= 2 print generated tables 


1. The SDATA card must be selected by the DLOAD Case Control card. 

2. This loading may be combined with other loadings by means of the DLOAD 
bulk data card. 

3. The items SAVE and SOURCE refer to the NASTRAN GINO file INPT, INP1 thru 
INP9. INPT is denoted by zero, INP1 - INP9 by the integers 1 to 9. 

4. If SAVE is blank, the output time histories will not be saved. 

5. If any DIRi or THNi is left blank, then the scan of these items is 
terminated. 

6. Up to 6 acceleration time histories at a single grid point may be 
defined on one logical card. 


87 



Input data card SET1 


Description: Defines a set of grid points at which output is desired. 

Format and example 


123456789 10 


SET1 

SID 

Gl 

G2 

G3 

G4 

G5 

G6 

G7 

+abc 


10 

6 









m 

19 

G8 

-etc- 



















Field Contents 

SID Set identification number (integer >- 0) 

Gl, G2 etc List of structural grid points (integer > 0 or "THRU") 

Remarks: 

1. These cards are referenced by the SDATA card. 

2. If "THRU" is used if must appear in field 4. Fields 6 to 10 will then 
be blank. 

3- All points referenced within a THRU list must exist. 
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GENERATE 

MATRICES 


Loop for 
more 

subcases 



Figure 1: Generation of Floor Spectra 

Problem Flow 
















ID A, B 

APP DMAP 
TIME 10 
BEGIN $ 

$INSERT DMAP SEQUENCE HERE 
END $ 

CEND 


Figure 5: Sample Executive Control Deck 


TITLE = 

TSTEP 

METHOD = 

SPC = 

DLOAD = 

FREQ = 

OUTPUT (XYPLOT) 
PLOTTER 


GENERATE FLOOR RESPONSE SPECTRA 
50 
10 
100 
200 
201 

NASTPLT MODEL D,0 


Figure 6: Sample Case Control Deck 
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BEGIN BULK 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

FREQ 

2 

.01 

.05 







FREQ 

201 

.1 

. 5 

r\ 

• ? 

1.2 

1 V 

J- • t 

2.0 

?.c; 

+ FI 

+F1 

3-0 

3-5 

4.0 

4.5 

5.0 


7.0 

8.0 

+F2 

+F2 

10.0 

12.0 

14.0 

20.0 

28.0 

30.0 

35.0 



SDATA 

200 

-1 

3 

2 

-1 


BRTEH 

2 

+SD1 

+SD1 

1 

2 








SET1 

1 

6 








TSTEP 

50 

3000 

.005 

2 






$ 










$ 

ALL OT 

HER GEOMETRIC DATA 







$ 










ENDDATA 











Figure 7: Sample Bulk Data Deck 
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zo mi ;d son r m o o 




Figure 9: Generated Floor Response Spectra - Printed 


GENERATION 

OF FLOOR 

RESPONSE SPECTRA 

SAMPLE PROBLEM 

JULY 

14 • 1380 

SINGLE POINT EXCITATION PROBLEM 



OUTPUT FOR 

COLLOQUIUM PAPER 





POINT-ID = 

6 

DIRECTION = 2 






F 

L 0 0 R R 

ESPON 

SE SPECTRA 



*• 



■ EQUIPMENT 

DAMPINGS 

— — — — - 

FREQUENCY 

TYPE 

•010 

-050 




3-200000+00 

G 

3-499337-01 

2-088033-01 

• 0 

• 0 

• 0 



• 0 

• 0 

• 0 

-0 

-0 

3-6D0D0D+00 

G 

3-340834-01 

1-794977-01 

-0 

• 0 

• 0 



• 0 

• 0 

• 0 

-0 

• 0 

3-800000+00 

G 

3-436410-01 

2-215128-01 

• 0 

-0 

• 0 



• 0 

• 0 

-0 

• 0 

• 0 

4 • 000000+00 

G 

3-850916-01 

2-188926-01 

• 0 

• 0 

-0 



• 0 

-0 

• 0 

• 0 

• 0 

4-200000+00 

G 

4-236383-01 

2-059928-01 

-0 

-0 

-0 



• 0 

-0 

-0 

-0 

-0 

4-800000+00 

G 

3-487423-01 

2-203638-01 

• 0 

• 0 

-0 



• 0 

• 0 

• 0 

-0 

-0 

5-200000+00 

G 

6-550089-01 

2-713436-01 

• 0 

• 0 

-0 



-0 

• 0 

• 0 

-0 

-0 

5-600000+00 

G 

5-732045-01 

2-750139-01 

• 0 

-0 

-0 



-0 

-0 

• 0 

• 0 

-0 

6-200000+00 

G 

1-040017+00 

4-595057-01 

• 0 

-0 

-0 



-0 

• 0 

• 0 

• 0 

-0 

6-600000+00 

G 

1-287466+00 

5-287637-01 

-0 

-0 

• 0 



• 0 

• 0 

• 0 

-0 

-0 

6-680144+00 

G 

9-027770-01 

5-140159-01 

• 0 

• 0 

• 0 



-0 

• 0 

-0 

-0 

• 0 

6-8 00000+00 

G 

8 - 883788-01 

4 - 994227-01 

• 0 

• 0 

-0 



-0 

•0 

• 0 

• 0 

• 0 

6-971636+00 

G 

3-545243-01 

4-633139-01 

• 0 

• 0 

• 0 



• 0 

-0 

-0 

• 0 

• 0 

7-000000+00 

G 

8-189829-01 

4-544786-01 

-0 

• 0 

• 0 



• 0 

• 0 

-0 

-0 

• 0 

7-600000+00 

G 

8-235375-01 

3-937230-01 

• 0 

• 0 

-0 



• 0 

• 0 

• 0 

-0 

• 0 
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SOLUTION OF ENFORCED BOUNDARY 
MOTION IN DIRECT TRANSIENT AND 
HARMONIC PROBLEMS 
Prepared By 
Gary L. Fox 
Director 

Engineering Mechanics Division 

NKF ENGINEERING ASSOCIATES, INC. 
8150 Leesburg Pike 
Vienna, Virginia 22180 
(703) 442-8900 

INTRODUCTION 


The current versions of NASTRAN, i.e., NASA, MSC, and MAC support non-zero 
boundary displacements only in the static analysis. Forcing functions in the 
dynamic analysis formats allow only forces and pressures to exercise the 
mathematical model. This limitation can be circumvented by the application 
of a DMAP alter sequence. For the direct harmonic problem, a simple change to 
module FRRD can be easily incorporated to effect a more efficient use of the 
code. 


Let the equation of motion be written with the dynamic set of coordinates 
in partition form with subscript b as the boundary set and subscript c as the 
complimentary boundary set, i.e.. 


where 


m m , 
cc cb 


^c 





cb 


bb 


1 

• \ 

r x J r 

c I 


• / + 

) 

-M 


cc 







m, d, k = mass, damping, and stiffness matrix coefficients 


P 2. = linear and non-linear load vectors 


96 



Equation (1) is not solved^by the direct transient or frequency formats when p 
X b , and therefore X^ and X^, are known and P^, X^, and therefore X £ and X^ c 
afe unknown. However, equation (1) can be rewritten in the form needed for 
solution by the standard NASTRAN modules. The first of these are: 

t" cc ] [“cel tv - tVMV (2) 

where 

t p c> - tP c > * [. cb ] [‘v * [d cb ] (X) ♦ [k cb ] [X b > 

By the use of the partitioning modules, the submatrices in Equations (1) or (2) 
are easily formed. By letting the boundary displacement vector be input 
through the FORCE or DLOAD cards, the force vector is actually identified as 
[P^} = [X,} (or the first or second derivatives). 

The formation of the load vector is different for the transient and har- 
monic cases. These issues will be discussed below. Somewhat independent of 
the problem is the requirement that the solution vector to be processed by the 
remaining modules must be of the dimensions of the "d" set. By using once 
more partitioning vectors and the MERGE module, the solution vector [X }, and 
in the transient case [ Y } and [X } , is merged with the boundary vecto? [X^l 
to form the dynamic vectSr [X^} . C With the "d" set solution vectors formed, 
the remaining DMAP sequence can be executed without NASTRAN knowing the 
difference. 

In the case of 
(2) becomes 

(w 1 2 [m cc ] + 

where 

w = circular frequency, 2 x f. 

HARMONIC ANALYSIS 


The DMAP alter that was written to partition the matrix equation (1) 
into the form of equation (2) and then solve the lower order equation (3) 
is shown in Figure A-l. The following paragraphs discuss the steps involved. 

1. FRRD calculates the load vector PDF and exits the module. The 
parameter ISKP is changed from -1 to a positive number to be 
transferred to FRRD the second time the module is executed. 

If the value of ISKP was set to zero, the default value, 
the module would have been executed normally. A normal 


harmonic analysis the non-linear force is zero and equation 
iw [d cc ] * [k cc ] ) [X c T = [P c l (3) 



execution would give a solution to equation (1). The FORTRAN 
listing of module FRRD is shown in Figure A-2. The added code 
is underlined. Only the subroutines FRRD1A and FRRD1B are 
executed in this step. 


2., The parameter ISKIP is saved for later use. 

3. The partition vector DPAR is used to partition the stiffness 
matrix KDD. The submatrix identification is related to equa- 
tion (2) by the following: 

Figure A-l. DMAP Alter for Harmonic Response 


ALTER 159.159 

FRRD CASEXX, USETD, DLT, FRL, GMD, GOD, KDD, BDD,MDD,DIT/UDVF,PSF, PDF, PPF/ (1) 

C, N, DISP/ C, N, DIRECT/ V,N,LUSETD/V,N,MPCF1/V,N, SINGLE/ V,N, OMIT/ 

V, N, NONCUP/ V, N, FRQSET/ V, N, ISKIP=-1/ $ 

SAVE ISKIP $ (2) 

PARTN KDD,DPAR,/KD11,KD21,KD12,KD22/ $ (3) 

PARTN MDD,DPAR,/MD11,MD21,MD12,MD22/ $ (4) 

PARTN PDF, ,DPAR/PD11,PD21,PD12,PD22/C,N,1 $ (5) 

MPYAD KD11 ,PD21,PD11/P1DF/C,N,0/ C ,N,-1/ $ (6) 

FRRD CASEXX, USETD , DLT , FRL , GMD , GOD , GOD , KD11 , , MD11 , DIT/UIDVF , PSF , P1DF , ( 7 ) 

PPF/ C,N,DISP/C,N, DIRECT/ V,N,LUSETD/V,N,MPCF1/V,N,MPCF1/V,N, SINGLE/ 

V, N, OMIT/ V,N, NONCUP/ V,N,FROSET/V,N, ISKIP/ $ 

MERGE KD11,KD21,KD12,KD22,DPAR,/KDD/ $ (8) 

MERGE MD11,MD21,MD12,MD22, DPAR, /MOD/ $ (9) 

MERGE U1DVF,PD21,PD22, ,OPAR/UDVF/C,N,l $ (10) 


ENDALTER 




CEND 






Figure A-2. Listing of Module FRRD 




LEVEL 2.2.1 (DFC 77) 


ISN 0002 


SUBROUTINE FRRD 

00000010 


C 


00000020 


C 

FREQUENCY AND RANDOM RESPONSE MODULE 

00000030 


C 


00000040 


C 

INPUTS CASECC, USETD, ULT, FRL, GMD, GOD, KDD, 




BCD , MDD , PHIDH , DIT 

00000050 


C 


00000060 


C 

OUTPUTS UDV,PS ,PD,MP 

00000070 


C 


00000080 


C 

8 SCRATCHES 

00000090 


C 


00000100 

ISN 0003 


INTEGER SINGLE , 0NIT , CASECC , USETD , DLT , FRL , 




GMD, GOD, BDD, PHIDH, DIT, 1 SCR1,SCR2,SCR3, 




SCR5 , SCR6 , UDV , PS , PD , FP , PDD , OPTION 

00000120 

ISN 0004 


INTEGER SCP7 , SCRB,NAME&2< ,MCB&7> 

00000130 



Oft 




ISN 

0006 


INTEGER FOL 

00000140 

ISN 

0006 


COMMON/ APP&2< ,M0DAL&2< ,LUSETD, MULTI .SINGLE, 





OMIT, NONCUP, FRQSET, 

00000150 



1 

ISKIP 

00000155 



C 


00000160 

ISN 

0007 


COMMON/FRRDST/ OVF&150< , ICNT , IFRST, ITL&3< IDIT, 





IFRD,K4DD 

00000170 

ISN 

0008 


DATA CASECC , USETD , DLT , FRL , GMD , GOD , KDD , HDD , 





MOD,PHIDH,DIT/ 

00000180 




1 101,102,103,104,105,106,107,108,109,110,111/ 

00000190 

ISN 

0009 


DATA UDV,PS.PD,PP,PDD/201,202,203,204,203/ 

00000200 

ISN 

0010 


DATA SCR1, SCH2 , SCR3, SCR4 , SCR5 , SCR6 /301,302, 





303,304,305,306/ 

00000210 

ISN 

0011 


DATA SCR7,SSCR8/307. 308 / 

00000220 

ISN 

0012 


DATA MODA /4HM0DA/ 

00000230 

ISN 

0013 


DATA POL/205/ 

00000240 

ISN 

0014 


DATA NAME /4HFRRD,4H / 

00000250 



C 


00000260 



C 

BUILD LOADS ON P SET ORDER IS ALL FREQ. 





FOR LOAD TOGETHER 

00000270 



C 


00000280, 

ISN 

0015 


IF ( ISKIP .GE. 0 ) GO TO 5 

00000281 

ISN 

0017 


NLOAD = ISKIP / 2*16 

00000282 

ISN 

0018 


NFREQ = ISKIP - NLOAD/ *2 **16 

00000283 

ISN 

0019 


GO TO 15 

00000284 

ISN 

0020 

5 

CONTINUE 

00000285 

ISN 

0021 


CALL FRPDIA&DIT , FRL , CASECC , DIT , PF , LUSETD , 





NFREQ , NLOAD , FRQSET , FOL , 

00000290 




1 NOTRD< 

00000300 

ISN 

0022 


1F&MULTI.LT.O. AND. SINGLE ,LT.O. AND. OMIT. L.T.O 

00000310 




AND. MODAL 





1 & 1< .NF. M0DA< GO TO 60 

00000320 



C 


00000330 



C 

REDUCE LOADS TO D OR H SET 

00000340 



C 


00000350 

ISN 

0024 


CALL FPRU14$PP .USETD , GMD , MULTI , SINGLE, OMIT , 

00000360 




M0DAL&1< , PH1DH , PD , 





1 PS , SCR5 , SCR1 , SCR2 , SCR3 , SCR4< 

00000370 

ISN 

0025 


IF ( ISKIP .LT. 0 ) GO TO 40 

00000375 

ISN 

0027 


15 CONTINUE 

00000377 

ISN 

0028 


IF ( MULTI .LT. 0 .AND. SINGLE. LT.O .AND. 





OMIT .LT. 0 

00000378 




. .AND. MODAL (1) .NE. MODA ) POD = PD 

00000379 



C 


00000380 



C 

SCR5 HAS PH IF MODAL FORMULATION 

00000390 



C 


00000400 

ISN 

0030 


IF &M0DAL&1< . EQ . MODA< PDD #SCR5 

00000410 



C 


00000420 



C 

SOLVE PROBLEM FOR EACH FREQUENCY 

00000430 



C 


00000440 

ISN 

0032 


IF&NONCUP .LT. 0 .AND. M0DAL&1< .EQ. MODA< 

00000450 


99 



GO TO 50 


ISN 

0034 


10 IF&FREQ .EQ. 1 .OR. NLOAD .EQ 1< SCR6 # UDV 

00000460 

ISN 

0036 


DO 20 1//1,NFPEQ 

00000470 

ISN 

0037 


CALL KLOCK&LOCK&ITIMEK 

00000480 



C 


00000490 



C 

FORM AND DECOMPOSE MATRICES 

00000500 



C 

' 

00000510 

ISN 

0038 


CALL FRRD1C&FRL,FR0SET,MDD,RDD,KDD.1,SCR1, 





SCR2 , SCR3 , SCR4 , SCR8 , 





1 SCP7.1G00D< 

00000530 



r 

Vj 


00000540 



c 

ULL IS ON SCR1 — LLL IS IN SCR2 

00000550 



c 


00000560 



c 

SOLVE FOR PD LOADS STACK ON SCR6 

00000570 



c 


00000580 



c 


00000590 

ISN 

0039 


CALL FRRD1D&PDD , SCR1 , SCR2 , SCR3 , SCR4 , SCR6 , 





NLOAD , 1G00D , NFREQ< 

00000600 

ISN 

0040 


CALL KL0CK&ITIME2< 

00000610 

ISN 

0041 . 


CALL IMTOGO&ITLEFK 

00060620 

ISN 

0042 . 


IF&2*&ITIME2-ITIME1< .GT. ITLEFT .AND. I .NE. 

00000630 




NFREQ< GO TO 70 


ISN 

0044 

20 

CONTINUE 

00000640 

ISN 

0045 


1 # NFREQ 

00000650 

ISN 

0046 

30 

CONTINUE 

00000660 

ISN 

0047 


IF&NFREQ .EQ. 1 .OR. NLOAD .EQ 1< GO TO 40 

00000670 



C 


00000680 



C 

RESORT SOLUTION VECTORS INTO SAME ORDER AS LOADS 

00000690 



C 


00000700 

ISN 

0049 


CALL FRRD1E&SCR6, UDV, NLOAD, I< 

00000710 

ISN 

0050 


40 ISKIP = NFREQ +NLOAD*2**16 

00000720 

ISN 

0051 


RETURN 

00000725 



C 


00000730 



C 

UNCOUPLED MODAL 

00000740 



C 


00000750 

ISN 

0052 


50 CALL FRRD1F&MDD, HDD, KDD,FRL,FRQSET, NLOAD, 





NFREQ, PDD,UDV< 

00000760 

ISN 

0053 


GO TO 40 

00000770 

ISN 

0054 


60 PDD # PP 

00000780 

ISN 

0055 


GO TO 10 

00000790 



C 


00000800 



C 

INSUFFICIENT TIME TO COMPLETE ANOTHER LOOP 

00000810 

ISN 

0056 


70 CALL MES AGE& . 5 . NFREQ-I , NAME< 

00000820 

ISN 

0057 


MCA&1< # SCR6 

00000830 

ISN 

0058 


CALL RDTFL&MC A*l< < 

00000840 

ISN 

0059 


MDONE # MCD&2< 

00000850 

ISN 

0060 


MCR&1< # PP 

00000860 

ISN 

0061 


CALL R0TR1&MCH&1« 

00000870 

ISN 

0062 


MCR&2< NOONF 

00000880 

ISN 

0063 


CALL WRT1FL&MCB&1 

00000890 



ISN 

ISN 

ISN 


0064 

IF&SINGLE .LT. 0< GO TO 80 

00000900 

0066 

MCA&1< # PS 

00000910 

0067 

CALL PUTRL&MCA&1« 

K, . = KD11 
K d * = KD12 
Kf b = KD21 

Sb " 022 

00000920 




4. The partition of the mass matrix, MDD, is similar to the stiffness 
matrix. 

5. Because the load vector is calculated for all frequencies and 
loading conditions at once, PDF is a load matrix, a load vector 
in each 'Column. The partition vector DPAR is used again to 
separate the enforced displacements from the forces. The rela- 
tionship to equation (2) is 

P = PD11 

P^ = PD21 
b 

6. The module MPYAD performs the matrix multiplication and additions 
required by equation (2). Here 

P = P1DF 
c 

7. Module FRED is executed again, but this time the parameter ISKIP is 
positive. A jump to statement 15, underlined in Figure A-2, causes 
only the subroutines FRED1C, FRRD1E and FRRD1F to be executed. The 
solution to equation (3) is obtained in this step. The code uses 
the following names related to equation (3). 

M = MD11 

K CC = KD11 

P C = P1DF 

X C = U1DF 
c 

8. The stiffness matrices are merged to form the system stiffness matrix. 
This is the inverse of operation 3. 

9. Similar to the stiffness matrix, this operation is the inverse of 
operation 4. 

10. Merges the solution vector X of equation (6-7) with X^ to form the 

system solution vector X,. 

d 

The three merges, operations 8, 9, and 101, are made necessary because 
NASTRAN uses the displacement approach to the problem solution. In 
order to calculate stress and forces in the members, the solution 
vector must contain all grid points. 
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TRANSIENT ANALYSIS 


The DMAP Alter required for the Rigid Format 9, Direct Transient Response, 
is shown in Figure A-3. The discussions below relates to the circled numbers 
in the DMAP listing. 

1. The Stiffness matrix is partitioned in accordance with 
Equation (2) where 

KD11 
KD12 
KD21 
KD22 

2. The Mass matrix is partitioned similar to the Stiffness matrix 

MD11 ! MD12 
i 

r — 

MD21 | MD22 

Figure A-3. DMAP Alter to Rigid Format 9 


ALTER 163 

PARTN KDD , OPAR, / KD11 , KD21 , KD12 , KD2 2 / $ (1) 

PARTN MDD,OPAR,/MDLL,MD21,MD12,MD22/ $ (2) 

PARTN PD,0PAR/PD11,PA21,PD12,PD22/C,N,1 $ (3) • 

MPYAD PA21,MV1,/PBT21/C,N,0/C,N,1/C,N,0/C,N,2 $ (4) (5) 

ADD PBT21 , PA21/PB21/C ,Y , ALPHA= (0. 550E-2 . 0) /C,Y,BETA= (0 . 550E-2 . 0) $ (6) 

MPYAD PB21,MAIT,/PV21/C,N,0/C,N,1/C,N.0/C,N,2 $ (7) 

MPYAD PV21,MV1,/PCT21/C,N,0/C,N,1/C,N,0/C,N,2 $ (8) 

ADD PCT21 ,PV21/PC21/C,Y, ALPHA= (0. 550E-2 . 0. ) /C,Y,BETA= (0. 550E-2 . 0) $ (9) 

MPYAD PC21,MAIT,/PU21/C,N,0/C,N,1/C,N,0/C,N,2 $ (10) 

MPYAD KD12,PU21,PD11/P1D/C,N,0/C,N,1 /$ (11) 

ALTER 165,165 

TRD CASEXX,TRL,NLFT,DIT,KD11,MD11,PID/UIDVT,P1LD/C,N, DIRECT/ 

V,N,N0UE/V,N, NONCUP/ V,N,NC0L $ (12) 

ALTER 166 

MERGE KD11,KD21,KD12,KD22,0PAR/KDD/ $ (13) 

MERGE MD11 ,MD21 ,MD12 ,MD22 , OPAR, /MDD/ $ (14) 

MERGE PD11 , PILD , PD12 , PD2 2 , , 0PAR/PNLD/ C , N , 1 $ (15) 

PARTN PA21,PVA,/A21,,PDA12,/C,N,1 $ (16) 

PARTN PV21,PVA,/V21,,PDA12,/C,N,1 $ (17) 

PARTN PU21,PVA,/U21,,PDA12,/C,N,1 $ (18) 

MERGE A21,,V21,,PWA,/PVA21/C,N,1 $ (19) 

MERGE PVA21, ,U21, ,PVUVA, /PUVA21/C,N,1 $ (20) 

MERGE U1DVT,PUVA21, , , ,DPAR/UDVT/C,N,1 $ (21) 

END ALTER 
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3. The load vector, PD, which is output from module TR LG, is partitioned 
according to Equation (2), where 

PD = (P(t 1 )>, {P(t 2 )}, ... 

PD11 = (P (t n )}, {P (t_)}, . . . 
cl c l 

PA21 = (P b (t 1 )}, (P b (t 2 )}, . . . 


Note that PD is a matrix formed by columns of load vectors, one 
column for each time step. The matrices PD22 and PD12 are not 
generated, i.e. 


PD 


PDli 

PA21J 


4. Direct input matrices, MV1 and MA1T, are used subsequently to calcu- 
late the velocity and displacement matrices from the acceleration 
matrix. The forms of MV1 AND MA1T are 


MV1 


MA1T 


oioo..: 

1 

0 0 10. . . 


0 0 0 1 . . . 

1 

• • • • 

• • • • 

M 

1 

• • • • 

1 


— N+2— - 


'mi... 


0 111... 


0 0 11... 
• • • • 

M 

• • • • 

• • • • 

•* — N+2 



The dimensions of both matrices are M X N + 2 where M is the number 
of coordinates in the b-set and N is the number of time steps. 

5. Produces the matrix product 

[PBT21] = [PA21]*[MV1] 

= [{P b (t b )}, {P b (t 2 )}, . . .] 


"oioo..: 
0 0 10 . 

0 0 0 1 ... 


= [0, (P b 0^)}, {P b (t 2 )}, . . .] 


It is seen that this operation moves the columns of the acceleration 



vectors from time t. to t.+l. 

11 

6. Produces the matrix sum 

[PB21] = a [PBT21] + 8 [PAZ1] 

The coefficients a and 3 are set equal to one-half of the 
integration time step, Atl 

[PB21] = [{P 1 + P 2 >, {P 2 + P 3 >, . . .] 

=[{AV 1 >, (AV 2 >, . . .] 

where [P J = { P (t . ) }: i = 1 to N + 2 
1 c 1 

The matrix PB21 represents the change in velocity, AV., between 
time steps, t^ and t^ The calculation is based on the trapesoidal 
rule for numerical integration. 

7. The final step in producing the matrix of velocity vectors, PV21 
from the matrix of acceleration vectors, PA21, this module produces 
the matrix product 

[PV21] = [PB21] [MA1T] 

fun. . 

0 111. . 

0 0 11. . 

= [AV } , (AV }, . . .] 0 0 0 1 .. 


= [{AV X >, (AV 1 + AV 2 >, {A x + A 2 + AV 3 >, . . .] 

8., 9. A repeat of operations e, f, g. The matrix of displacement 
and 10. PU21, is calculated from the matrix of velocity vectors, PV21. 

11. The load vector is calculated in accordance with Equation (2). 

KD12 = K 

cb 

PU21 = {X b > 1 , {X b ) 2 ,. . . . 

PD11 = {P^, {P b > 2 ’ * * * 

P1D -{P}.,{P} 0 , ... 
cl c 2 

12. The module 'TRD calculates the solution to Equation (2). 

KD11 = [Kcc] 

MD1 1 = [Mcc] 

P1D = P , P 2 , . . . 
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U1DVT = 


X X 

• • 

X X 

• • • • 

X J, X 2 , 

[PUD] = {P^J . 

The solution vector, U1DVT, is a matrix of displacements, velocity 
and acceleration vectors for each grid point; a column for each 
time step. 



13. The system stiffness matrix is formed 


KD11 KD12 

KD21 KD22 


= [KDD] 


14. The system mass matrix is formed similar to the operation (13.) 


15. The system load vector is formed 


16., 17. 
and 18. 


PD11 

P1LD 


[PNLD] 


Partition the acceleration, PA21, velocity, PV21, and dis- 
placement, PU21, matrices to the correct size to be merged 
with U1DVT. 


19., 20. These operations merge the solution matrix, UD1VT, with the 
and 21. excitation matrix, PUVA21, to form the final system solution 
matrix, UDVT. 


UD1VT 


PUVA21 


[UDVT] 


From this point on, the solution is calculated according 
to the Standard Rigid Format 9 procedure. 
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APPLICATION OF NASTRAN TO THERMAL TRANSIENT ANALYSIS 
WITH SURFACE ABLATION* 

Karl Meyer 

Planning Research Corporation 
John F. Kennedy Space Center, Florida 

SUMMARY 


This paper presents a computer modeling technique for solving thermal 
transient analysis {Solution 9, Approach Heat) with surface ablation problems 
using the NASTRAN Computer Program. 

The mathematical model used in this analysis is one dimensional, which 
corresponds to the direction of heat flow. All dimensions perpendicular to 
that of the heat flow direction are assumed to be in thermal equilibrium, 
i.e., the net heat flow in and out is assumed to be zero. A special modeling 
technique that would simulate the effects of surface ablation corresponding to 
the direction of heat flow was developed and incorporated into this mathemat- 
ical model. 

All computer analyses were done using Level 16.0 NASTRAN on a UNIVAC Sys- 
tem. This analysis technique was developed to predict the thermal response 
and the amount of surface ablation that would occur on the deck of the Mobile 
Launcher Platform (MLP) during a launch of the Space Shuttle at John F. 
Kennedy Space Center, Florida. 


INTRODUCTION 


The similarities between heat transfer and structural analysis have been 
exploited to include heat transfer as a part of NASTRAN analysis capabilities. 

In heat transfer, just as in structural mechanics, the analysis of a 
solid continuum can be modeled using the finite element method. Finite ele- 
ment analysis reduces the problem to the solution of a set of simultaneous 
equations in which the unknown variables to be calculated are defined at a 
discrete or finite set of grid points. This set of simultaneous equations may 
then be expressed in the matrix form: 


[K] {u} + [B] {u} = {P} + {N} (ref. 1) 


*This work was done under NASA contract No. NAS10-8525. 
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Table I shows the analogies between heat transfer and structural mechan- 
ics. The major difference between heat transfer and structural mechanics 
analysis is that temperature is a scalar function of position, whereas dis- 
placement is a vector with as many as six component degrees of freedom. Thus, 
in heat transfer analysis, NASTRAN provides one degree of freedom at each grid 
point (temperature). 

The heat conduction matrix [K], and the heat capacity matrix [B], are 
derived from element and material properties assigned by the programmer. An 
additional advantage of NASTRAN is that part of the heat conduction matrix may 
be associated with surface heat convection or radiation. The applied heat 
flux vectors {P} , are associated with heat flux in or out and are time- 
dependent loads. The nonlinear heat flux vectors (N) are temperature depend- 
ent or rate-of-change-of-temperature dependent heat flux in or out of the sys- 
tem. 


Typical output from a NASTRAN thermal transient analysis includes: 

a. Temperature vs time or grid points 

b. Rate of change of temperature vs time at grid points 

c. Thermal flux vs time at grid points 

d. Finite element temperature gradient vs time 


[B] 

[K] 

(NT 

{P} 

{u} 

(u> 

A 

F 

S 


SYMBOLS 

Matrix of constant heat capacitance coefficients 

Matrix of constant heat conductance coefficients 

Matrix of nonlinear thermal flux which are functions of tempera- 
ture of the rate of change of temperature 

Matrix of thermal flux are functions of time. 

Matrix of temperature at grid points 

Matrix of the rate of change of temperatures with respect to time 
at grid points 

Unit cross-sectional area 

Radiation factor 

Arbitrary scale factor 
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b 

c 

h 

k 

AX 

p- 

<a 

f( ) 

c 

i 

in 

out 

TS 

1 , 2 , etc. 


Thermal capacity 
Specific heat 

Convective film coefficient 
Thermal conductivity 

Length of element x = .0529 cm (1/48 inch) 
Density 

Stefan-Boltzmann constant 

Function of a variable 

Second component CELASi element or point 

The i^h grid point or element 

Thermal flux in 

Thermal flux out 

Thermal switch 

Grid point or finite element number 


NASTRAN MODELING PROBLEM 


When analyzing a thermal transient heat transfer problem using a finite 
element program such as NASTRAN' in which surface ablation is present, a 
modeling problem occurs. The problem is that when an element is ablated from 
the material being analyzed it must be eliminated from the analysis. In the 
present configuration of NASTRAN, it is not possible to eliminate an element 
or elements in the middle of the computations and continue on to the 
completion of the computations for a complete thermal transient analysis. 

To solve this problem, a NASTRAN computer model was developed consisting 
of a one-dimensional isotropic material or materials in which transient 
transport of thermal energy occurs, with ablation from the front surface. All 
dimensions perpendicular to that of heat flow direction consist of isotropic 
material (s) which is assumed to be in thermal equilibrium, i.e., the net heat 
flow in and out is assumed to be zero. 

This NASTRAN model uses a special modeling 'technique to simulate the 
removal of elements caused by surface ablation. The modeling technique 
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involves the use of NASTRAN nonlinear load bulk data cards. These nonlinear 
loads essentially cancel an element from the analysis when it is ablated and 
relocate the thermal flux into the model to the next prescribed grid point. 


NASTRAN THERMAL MODEL 


Refer to figure 1 (NASTRAN Thermal Model) in the following discussion. 


To model the thermal properties of the material being analyzed, NASTRAN 
CROD. elements connected in series were used. This series connection of CROD 
elements describes the one-dimensional thermal properties of the material or 
materials being analyzed with respect to depth. For most heat conduction 
problems, 20 elements per centimeter (48 elements per inch) of thickness are 
needed to represent accurately the component being analyzed (ref. 2). Refer- 
ing to figure 1, the thermal conductors Kj to K(-j_j) and thermal capacitors 
Bj to B-j are modeled using CROD elements 1 to (i - 1). The CROD elements ref- 
erence PROD property cards in which a unit area (A-j) of one was used. The 
PROD property card references MAT4 material property card. The MAT4 material 
property card requires the thermal conductivity (k-j) and the thermal capacity 
(b-j) of the material (s) being analyzed. NASTRAN calculates the thermal con- 
ductance and thermal capacitance for each element using the following rela- 
tionships: 

ftp = k-j 'A-j/ AX 


and B-j = bjA-jAx/2 (lumped method) 


where bf = pic-j 


The backface grid point -j has a CHBDY element (convective film coeffi- 
cient h) which references a temperature at grid point r. The CHBDY element 
references a PHBDY property card where an area factor of one is specified. 
The PHBDY card references a MAT4 material where the convective film coeffi- 
cient is specified. 

To model the ablation of elements from the model, a concept was developed 
which I will call a thermal switch. The thermal switch is an element that has 
an output response of zero or one, depending upon some predefined ablation 



temperature at a reference grid point in the material model. The thermal 
switch was modeled using grounded CELAS3 elements. There is one thermal 
switch for each point in material model. These thermal switches are not 
physically connected to the grid points in the material model. The response 
temperature of each thermal switch is made a function of the temperature at a 
reference grid point using N0LIN1 nonlinear loads. The N0LIN1 nonlinear load 
applies a thermal load N at the ungrounded point of the CELAS3 element, which 
is a function of a temperature at the reference grid point and a table 
(TABLEDi) bulk data card. The response temperature at the ungrounded point of 
the CELAS3 is given by the following relationship: 



It is evident from this equation that if a thermal load N of zero is 
applied when the temperature at the reference grid point is below the ablation 
temperature, the temperature response u will be zero. When the temperature of 
the reference grid point is equal to or greater than the ablation temperature, 
a thermal load of N = K is applied. The resulting response temperature u at 
the ungrounded point of the CELAS3 element will be equal to one. Essentially, 
this modeling technique makes the CELAS3 element an on/off switch which is a 
function of temperature at a reference grid point in the material being 
analyzed. 

When the response temperature of the CELAS3 element ungrounded point is 
used in conjunction with N0LIN2 nonlinear loads and other CELASi elements, all 
effects of ablation can be modeled. The response of the other CELASi 
ungrounded points are used as the second components, either temperature 
dependent or rate-of-change-of-temperature dependent, in the N0LIN2 nonlinear 
load cards. The first component response for N0LIN2 nonlinear load card is 
the response temperature of the ungrounded point of the respective thermal 
switch. The general equation for the N0LIN2 nonlinear load card is: 


N-j = S-j •f(uj$ i ).f(u cl - ) (ref. 3) 

where: S-j is an arbitrary scale factor for the jth point, f( ujs^ ) 
is the response of the jth thermal switch, f( u c j ) is the second 
component response of the jth CELASi element ungrounded point. 



The second component response of the unground point of the CELASi element 
can be used to generate thermal loads such as: 

f(u C i)=M ■ 

where, N c i = j^Fcru^ for radiation thermal loads in or out 
N c i = _+hu for convective thermal loads in or out 
N c i = ±Ku for conduction thermal loads in or out 
N c i = j^Bu for capacitance thermal loads in or out 


When the second component response is used in conjunction with the 
thermal switch response along with the N0LIN2 nonlinear load card, the 
resulting loads can be turned on or off based on some predefined ablation 
temperature. 

The ablation effects that can be modeled using this modeling technique 
include: 


a. Cancelling of the original thermal flux in and relocating the thermal 
flux to the next prescribed grid point 

b. Cancelling the thermal conductance effects of the ablated element by 
applying an equal and opposite thermal load at the next grid point 
of: 


Nout-j = -K i*( u i-1 “ u i ) ( u TS-j ) 

c. Cancelling the thermal capacitance effects of the ablated element by 
applying an equal and opposite thermal load at the next grid point 
of: 


^outf ~ ^i-l ( u TSj_2 ) 
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NOTE 


This effect is negligible and can be assumed to be 
zero, because the value of (ij-l approximately equals 
zero when the other ablation effects are taken into 
consideration. Also, this type of nonlinear load can- 
not be used in NASTRAN's present confi guration. 
(Refer to NASTRAN PROGRAM RECOMMENDATIONS section of 
this report.) 


Again, all these ablation effects can be modeled with thermal switches in 
combination with N0LIN1 through N0LIN4 nonlinear loads and other grounded 
CELASi elements. One should note that when CELASi elements are used to pro- 
vide non-zero temperature responses, a large thermal conductive to ground is 
required. For the relationship N = Ku, u must be adjusted to the desired tem- 
perature by defining the thermal conductance, K, of the CELASi element, which 
is connected to ground, and a load N, which is applied to the grid or scalar 
point in question. The numerical value of K should be several orders of mag- 
nitude greater than the numerical value of the thermal conductance specified 
in the rest of the model. 

Figure 2 shows a schematic diagram of all thermal load equations applied 
in a typical ablation program. The equations on the right of figure 2 are 
used for applying and cancelling thermal flux in after the ablation process 
occurs. Initially, N-j nl is applied at grid point 1. When grid point 1 reaches 

ablation temperature, N out ] is turned on by thermal switch 1 (CELAS3 element) 
and this thermal load is applied to grid point 1. The value of N outl is equal 
and opposite to that of N-j n -|, essentially cancelling the thermal flux into 
grid point 1. At the same time N ou f| is turned on Nj n 2 , which is equal to 
Ninl» is turned on, thus moving the input thermal flux from grid point 1 to 
grid point 2. This process is repeated for grid point 2 to grid point i, when 
and if each respective grid point reaches ablation temperature. 

The equation on the left of figure 2 are used for cancellation of element 
thermal effects after ablation of the respective element(s) occurs. Initi- 
ally, none of these thermal loads are applied until the ablation process 
occurs. When grid point 1 reaches ablation temperature, N ou ^2 is turned on by 
thermal switch 1 and this thermal load is applied to grid point 2. The N ou t2 
thermal flux cancels the thermal flux of Bj and Kj. This process is repeated 
for grid points 3 to grid point i, when and if each respective grid point 
reaches ablation temperature. 
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Multilayer surfaces are modeled in the same manner and are composed of 
layers of materials having different thicknesses and thermal properties. 


NASTRAN PROGRAM EXECUTION 


It was found during program execution that using a value of 6 = 1 for the 
Newmark Beta integration algorithm (ref. 1) proved to be the most efficient 
with respect to computer processing time. The value of ,$ is specified on a 
NASTRAN PARAM bulk data card. It should also be noted that the DIAG 10 option 
can be used in the executive control deck, which will result in the use of the 
alternate nonlinear load technique. This option may allow the user to use 
larger integration steps depending upon the particular problem. 


NASTRAN RESULTS 


Figures 3 and 4 show the theoretical thermal response plots of the 
NASTRAN thermal program using this modeling technique. These curves show the 
thermal response of the surface (Tsurface) exposed to the thermal flux in, at 
25% ( T 25%) » at 50 % ( T 50%) » at 75 * ( t 75%) and at the backface (T BACKFACE ) of a 
1-inch-thick steel plate. The flat portion at the top of the curves shows 
when ablation occurred at each point. The plots were generated using the 
NASTRAN X-Y plot capability. 

These two thermal response plots were used to verify the NASTRAN program 
against the actual physical results of two test situations. The tests con- 
sisted of the impingement of the plume of a Tomahawk missile normal to a 
2.54-centimeter (1-inch) thick steel plate for 9 seconds at two different 
relative positions. For the first test, the nozzle exit plane of the Tomahawk 
missile was located 5.182 meters (17 feet) from the steel plate. The results 
of this test showed 75% ablation occurred. This test corresponds to figure 3, 
which shows the same results. 

In the second test, the nozzle exit plane of the Tomahawk missile was 
located 3.658 meters (12 feet) from the steel plate. The results of this test 
shows a burn-through (100% ablation). This test corresponds to figure 4, 
which shows the same results. 

Figure 5 is included to show a thermal response comparison between run- 
ning NASTRAN without this modeling technique and with the modeling technique 
(figure 4). Large differences between figures 4 and 5 are quite apparent. 
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CONCLUSIONS 


The inclusion of nonlinear loads in the NASTRAN computer program allows 
the user a tremendous versatility in the solution of dynamic problems. When 
using these nonlinear loads in conjunction with NASTRAN, the user has a fast 
and comparatively easy tool in solving problems of thermal transient analysis 
with surface ablation. 

It should also be noted that this type of modeling technique could be 
applied to other types of dynamic problems when using NASTRAN. 


NASTRAN PROGRAM RECOMMENDATIONS 


For the solution of thermal transient response problems, as well as for 
other types of dynamic response problems, the following NASTRAN program recom- 
mendations are made: 

a. Acceleration-dependent nonlinear loads should be included in the 
NASTRAN NOLINi bulk data cards. 

b. All combinations of the three different dependencies, displacement, 
velocity, and acceleration, should be included in the NASTRAN N0LIN2 
nonlinear load bulk data card. 

The implementation of these recommendations into NASTRAN in conjunction 
with the on/off switch and existing NOLINi nonlinear load dependencies, would 
allow the user to simulate disconnections (uncoupling of the equations) and/or 
alter a section or sections of a problem during program computations and con- 
tinue on to the completion of computations for a complete transient or modal 
transient analysis. 
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TABLE I. NASTRAN ANALOGIES 


Variable 

Heater Transfer 
Representation 

Structural Mechanics 
Representation 

u 

Temperature 

Displacement 

• 

u 

Rate of change of 
temperature (du/dt) 

Velocity 

P,N 

Heat Flux 

Forces, Moments 

K 

Thermal, conductance 

Stiffness 

B 

Thermal capacitance 

Damping 
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F I GURE I . NA ST RAN THERMAL MODEL 




FIGURE 2. NASTRAN ABLATION EQUATIONS MODEL 


£QUAT!ONS£OR CANCELLATION 
OF ELEMENT THERMAL EFFECTS 
AFTER ABLATION OCCURS 


EQUATIQNSEOR-ARRLYING, 
AND CANCELING THERMAL 
FLUX IN AFTER ABLATION 
OCCURS 


1 1 1 

hJa 

j SURFACE EXPOSED TO 
f~ \ THERMAL FLUX IN 

N inf h; «u in ) 

1 ( rn i n pn i mt l 

< 

< 

B 2, 

> K i 

> 

i 

' N outj = ’ s l" ,(u in ) ‘ f<u TS 1 ,! 
1 Nin.-Sj- «U in )-f(U TSl ) | 

c f PPin POINT 9 i 

N ou* , ■ K ,u r u 2» - f l 
B l u |.)* f ,u TS !> 4 

>Ko 
> L 

\ oK 1 U rUliMI L ! 

Nout 2 “ ’Sl • «U in > * «U T s 2 >l 

1 N in 3 - S r f,U in>-« u TS 2 > ! 

i rnin po i mt ^ 

N out, = [ 'K 2 (U 2 -U 3 )- t ■ 

b 2 0 2 ] • nu TS 2 ) Bj ' 


f "out 3 --Sr « u in>-«U T S 3 » 
N, n i - 1 “ s r «U in )-f(U TS ._ 2 > ' 

... / PPin POINT I .1 

N out j.j 55 [~ K j_2 ^ U |-2" 

-B: A if(U T , ) < 
” 2 ,_2J TS i-2 Bii 

i-l 

L 

1 uKIU rUIIMI I“1 

^out i-l* 

|N in .“V f(Uj n ) • f(U TS j j) 

1 [ HR ID POINT i 

N , - [-K. (U. ,-U.i- ♦ ; 
out j L 1 i-l 1 

B. A f(U TC ) 1 

i-l i-l J TSj.j 


— " j uf\ iu rv/iiMi i 

N out j ■ -Sl * f(Uj n ) • fWJfs ,) 
BACKFACE 


117 



i FIGURE 3. NASTRAN X-Y PLOT (TOMAHAWK TEST#1 SIMULATION) 

MODELING TECHNIQUE FOR ABLATION USED ! 




FIGURE 4. NASTRA N X-Y PLO T (TOMAHAWK TEST #2 SIMULATION) 
MODELING TECHNIQUE FOR ABLATION USED 
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ON THE APPLICATION OF NONLINEAR LOAD ELEMENTS 


TO THERMAL ANALYSES USING THE NASTRAN THERMAL ANALYZER 

Hwa-Ping Lee 

NASA/ Goddard Space Flight Center 


SUMMARY 


Using the nonlinear load elements in thermal analysis to simulate an 
undocumented nonlinear thermal boundary condition is presented. The treatment 
of the nonlinearity arising from the temperature-dependent convective film 
coefficients is shown in detail. As an illustration, emphasis is placed on 
the modeling techniques and their interrelationships with the solution accu- 
racy as affected by a specific integration algorithm of the transient thermal 
analysis used in the NASTRAN Thermal Analyzer. Briefly shown is the under- 
lying theory on which the maneuvering of terms pertinent to the modeling 
depends. This demonstration provides some insight into the intricacies of the 
method that would be general to all applications. A recommendation is also 
made to modify a nonlinear load element that will enhance the solution capa- 
bility and broaden the scope of application. ' 


INTRODUCTION 


The known solution capabilities of the NASTRAN Thermal Analyzer (NTA) 
have been well documented (ref. 1 and 2), and a large part of them have been 
demonstrated in detail (ref. 3). Despite their frequent encounterment in 
engineering applications, certain nonlinear solution capabilities have not 
been provided by the NTA in the transient-state thermal analysis (APP HEAT, 

SOL 9). Those excluded are, for instance, the temperature-dependent thermo- 
physical properties, temperature-dependent convective film coefficients, etc. 
Their absence from the original NTA's capability list was solely attributed to 
the fact that the required computer processing time (CPU and 1-0 times) would 
be excessive. Computations involving the reassemblages of the thermal conduct- 
ance matrix, together with the decompositions of the relevant matrices that 
contain the nonlinear thermo physical properties, would have to be repeated for 
all steps of the time increment in integration. The program shortcomings, 
however, may be compensated by an appropriate application of the available 
nonlinear load elements (NLLEs) in NASTRAN, although some elements were origi- 
nally developed for structural analyses. The intrinsic modular structure of 
NASTRAN permits them to be accessible to the NTA. Those NLLEs, obviously, can 
function as nonlinear thermal loads, and they can be employed to simulate other 
complex physical phenomena as well. Otherwise, new modules would have to be 
developed and implemented, even though an excessive computer processing time 
was tolerable. Such an alternative would be an expensive and time-consuming 
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proposition. The user would be deprived of immediate solution meanwhile. 

Six NLLEs are currently available in NASTRAN, namely N0LIN1*, N0LIN2, 
N0LIN3, N0LIN4, N0LIN5 and NFTUBE. The last two were developed specifically 
for thermal analyses: N0LIN5 is capable of simulating temperature sensitive 

thermal radiation surface properties, such as a louver (ref. 4), and NFTUBE is 
capable of dealing with thermal energy transfer between the wall and the fluid 
flowing inside a tube (ref. 5). The remainders are four nonlinear elements of 
general type capable of accommodating four different functions with a temper- 
ature or heat flux being the argument. N0LIN1 allows to specify a table through 
TABLED! (i = 1,2, 3, 4), and it is the most versatile NLLE in NASTRAN. 

The main purpose of this paper is to illustrate how to acquire a non- 
routine, but essential solution capability by an application of the NLLEs 
through modeling. For a systematic demonstration, N0LIN1 was employed to 
simulate the temperature-dependent convective film coefficients, h(T). The 
importance of this problem is evident in view of a variety of applications 
ranging from a domestic heating system, a hot gas flowing through a nozzle, 
a fluid flow in a nuclear reactor, and a flow in a wind tunnel during the 
start-up to a thermal ablation including a receding surface. They are all'" 
characterized by such a nonlinear boundary condition. More specifically, 
treatments of the temperature-dependent film coefficients based on one and 
two reference temperatures are given in detail. The dilemma arising from the 
presence of the MPC to establish an average temperature of two reference 
temperatures and the reference with that average temperature by the NLLEs was 
resolved by a trick modeling. A brief review of the theoretical aspect serves 
not only to provide a basis for the maneuvering in modeling and to provide 
insight into the intricacies of the method, but to amplify basis of modeling 
techniques that are equally valid to other applications. A recommendation for 
the improvement of the N0LIN2 is outlined. It would enhance solution capa- 
bilities and broaden the scope of engineering application. Guidelines for the 
user to apply the NLLEs in thermal analysis are provided. 


THE DEMONSTRATIVE PROBLEM 


For the convenience of easy reference and discussion, and later on, to 
demonstrate the modeling techniques, a relatively simple problem without losing 
any generality has been designed, and it is described below: 

The transient temperature responses in a composite slab (pipe) of an infinite 
extent are to be determined. Initially, the slab has a uniform temperature of 
293 K (20 C) throughout its body. The boundary surface of the metallic layer 
is in contact with a stream of hot gas maintained at 1422°K (1149°C) from time 
t = 0. The dimensions and thermophysical properties of the slab are tabulated 
in Table 1. The system depiction, together with its finite element model, 
is shown in Figure 1. 


*The names of actual NTA cards are capitalized. 
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THE THEORETICAL CONSIDERATION 


The success of attempting a solution utilizing any NLLEs to simulate a problem 
with nonlinear thermophysical properties relies upon the following factors: 

(1) The nonlinear terms, arising from the temperature-dependent thermophysical 
properties, which are embraced in the square matrices of the thermal conduct- 
ance and/or the thermal capacitance, must be extractable and separable from 
their normal positions. The resulting products of the temperature-dependent 
icoef f icients (kA/£ an d hA or C) and their respectively associated variables 

i(T or t) must be transposed and merged into the nonlinear thermal load . 

vector . 

(2) An NLLE of an appropriate form suited to accommodate the transposed 
quantities must be selected. Any restraints associated with the use of the 
NLLEs must be strictly observed. 

(3) A full acquaintance with the characteristics of each NLLE and the 
integration algorithm for the transient heat transfer analysis (APP HEAT, SOL 9) 
is essential in controlling and accomplishing desired solution accuracy, 
stability and efficiency at wish. 

To serve as the theoretical basis for such a maneuvering, a brief review 
of pertinent mathematical operations is in order. 

A transient thermal analysis in the NTA is based on the general heat equation 
in the matrix form as follows (ref. 1): 

[ C ] { T } +[K]{T} = {Q f } + {Q n } (1) 

where 

T = the vector of temperatures at grid points 
T = the rate change of temperature vector 
K = the thermal conductance matrix of constant elements 
C = the thermal capacitance matrix of constant elements 
= the vector of applied linear thermal loads 
Q n = the vector of applied nonlinear thermal loads 

The convective coupling terms are normally incorporated with the thermal 
conductive couplings in the thermal conductance matrix [K] regardless of their 
dependency. When the convective film coefficient is temperature-dependent, the 
convective coupling terms embracing h(T) have to be separated from their orig- 
inal [ K„ J , and thus[K„] |Tj. The extracted terms being a new vector |h(T)| 
are transposed to the other side of the equal sign in equation (1) , and they 
become or merge with the vector of nonlinear thermal loads, {Q n }* The 
rearranged new equation conforms exactly to equation (1), meaning that it is 
then within the solution capability of the NTA. , Mathematically, it can be 
expressed by 
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( 2 ) 


[K 0 (k,h(T))][T} = [ K-l (k) ] {t } +{h (T)| 

The integration algorithm for equation (1) is based on the modified 
Newmark Beta Method (ref. 6). It is a forward and backward differencing 
method with a parameter 0 that enables the user to select a value in the 
range of O<0<1. When equation (1) is transformed into such a difference 
expression, it reads (ref. 1): 

(K]{dT„+l + <l-f»T ll }+ JJC]{t „ +1 - T n } 

- {uq ' +1 + a-tn {qA>} + a+» K}-f {Cl} «> 


the subscript n refers to the nth time step. The accuracy and stability of 
solution are closely related to 0 . Equation (3) can be rearranged to give 

[l(C] + fi[K]]{T n+1 }=[i[C] - (l-0)[K]]{T. n | 

At L At 

+ {eQn+l + (l-0)Qn} + a+0){Q£} -0{Qn-l} < 4 ) 

The mat ri x [1 / At C + 0 K] is d e c om posed i nto its t riangular^ factors^ and 
solved at each time step via the forward and backward substitution 
process. If [K] contains terms of k(T) or h(T), a search from the tables 
and replacement of new values based on {T n } have to be carried out for all 
terms residing inside those two square matrices that are associated with 
|T n +jJ and {T n |. Consequently, the matrix decomposition has to be performed 
for each time step. This requires a very excessive computer processing time. 
However, [K] needs be evaluated only once if its constituent terms are all 
constant. This situation is equally valid in [C] , although C(T) is not 
considered in this study. 

Relating the terms of nonlinear thermal load to the specific demonstrative 
problem, it is obtained from equations (1) and (2) that 


(5) 


( 6 ) 


{q 11 } = - [h(T>] 


J \ 


where 


0 
0 

' H 200' 


{ 


H x = h(T r )AT L 

h 200 = “ h( T r) AT 200 


124 



and h(T r ) denotes the dependency of the convective film coefficient in a refer- 
ence temperature T r which can be the wall temperature or an average temperature 
of the free stream and the wall temperature as will be discussed and treated in 
following sections. In the NTA, the N0LIN1 card defines an NLLE of the form 

q"(t) = s^Cr-jO:)) (7) 

where Q^(t) is the load being applied to the GRID i at time t, S-^ is a scale 
factor. F(Tj (t) denotes a tabulated function in Tj to accommodate H^(Tj) 
and ^(Tj) as defined in equation (6). The equivalent nonlinear thermal load 
functions are entered through a TABLEDi (i = 1,2, 3, 4) card. T (t) is any 
permissible time-varying temperature on which the functions depend. Other 
NLLEs, N0LIN2, N0LIN3 and N0LIN4, can be selected as long as their specific 
forms are suitable for respective applications. 


THE CONVECTIVE FILM COEFFICIENT h(T r ) 


The convective film coefficient h governs the magnitude of heat transfer 
between the hot gas and the boundary surface of structure, h, in general, can 
be expressed as h(T r ) where the reference temperature T r may refer to only one 
temperature, such as the wall temperature, or more than one temperature through 
an intermediary, such as the average value of the free stream temperature and 
the surface temperature of the wall. A typical example is a fully developed 
pipe flow. The expression for h is (ref. 7) 

h = 0.023 hz- (Re)°' 8 (Pr)°' 4 (8) 

Dh 


For large temperature differences, a mean temperature is used in evaluating 
all physical properties except the flow velocity. This procedure would involve 
very tedious interrations. To simplify computations, it is desirable to 
separate the temperature-dependent from the independent factors. An alterna- 
tive form is, therefore, obtained by an application of the equation of state 
and the continuity equation, which leads to 


h 


= 0.023 C r 


0. 2 - 0.6 0.8 - 0 . 2 . 0.8 

./!» Pr« G» Dh 0F> 

J- oo 


(9) 


As the free stream temperature T B is known, all quantities in the preceding 
equation except T w are constants, h is, therefore, a function of T w . For the 
problem of concern, it is h = h(T^). The values of this function obtained by 
using equation (9) are listed in Table 2. Similar expression, but referencing 
with an average value of the free stream temperature and the surface tempera- 
ture is also widely used in engineering practices. For this demonstrative 
problem, it is h = h(T r ) and T r = %(Tj + T 200)- The values of h vs T r are 
given in Table 3. 
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THE MODELING TECHNIQUES AND DISCUSSION 


A one-dimensional finite element model is sufficient to represent the 
demonstrative problem, and this model is shown in Figure 1, A listing of the 
input data deck is reproduced in Figure 2. With comments appropriately inserted, 
functions regarding various segments of the deck are self-explanatory. The 
space bounded by two rows of asterisks is reserved for adding a packet of cards 
which simulate the effect of the convective film coefficients. Three cases of 
h are considered: (1) h of a constant value, (2) h = h(Ti) , and (3) h =h(T r ) 

with T r = h(Ti + T200)- 

h of a constant Value 


To assure correctness of modeling and to evaluate the obtainable accuracy 
of solution by the nonlinear load approach, the case of a constant film co- 
efficient with h = 0.467 W/cm^-K was tested first. Because, this case is well 
within the documented capability of the NTA. It can be accomplished using a 
train of cards CHBDY, PHBDY and MAT4 by the conventional modeling method. The 
temperature results so obtained serve as a base to compare with that obtained 
by the nonlinear load approach. The two packets of input data cards, whose 
images are shown in Figures 3(a) and 3(b), were used in respective computer 
runs. The equivalent nonlinear thermal loads for the pair of NLLEs through 
N0LIN1 405 with TABLED1 4051 and TABLED1 4052 are given in Table 4. The effect 
of selecting the time-step size At for integration on solution accuracy was 
also studied. 

Temperature results summarized in Table 5 with the headings of Cases 1 and 
3 are solutions yielded by these two approaches. In both cases, the same time- 
step size (the At-set (a) in Figure 4) was used. The largest difference between 
the two temperature profiles occurs at t = 0.25 sec., suggesting that a refine- 
ment of the time-step size is warranted. To maximize the efficiency of the 
computer processing time, refinements of At were made only in the beginning 
segments of time. Three sets of time-step size, which are presented in the 
NTA TSTEP card format and identified as the At-sets (a), (b) and (c) , are shown 
in Figure 4. A total of five cases were computed, two (Cases 1 and 2) by the 
conventional modeling method and the other three (Cases 3, 10 and 11) by the 
nonlinear load approach. Temperature results are tabulated in Table 5 for 
comparison. It is evident that both modeling methods can yield comparable 
temperature solutions with a desired accuracy, provided an appropriate At-set 
is used (e.g. , to compare Case 2 with Case 11). 

This test gave confidence in the nonlinear load approach and verified the 
adequacy of the two sets of time-step size, namely the At-sets (b) and (c) . 

They were, therefore, selected in all following computations. Inevitably, 
efficiencies of those runs using the At-sets (b) and (C) were penalized by a 
certain amount. They are reflected by the CPU and 1-0 times, which are also 
summarized in Table 5. 
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h = h(TT) 


The modeling procedure similar to that used in the preceding case of a 
constant h is equally applicable to this case of a convective film coefficient 
h, which is dependent on the surface temperature of the wall. The preparation 
of the equivalent nonlinear thermal loads, however, is much more involved. 
Computations of those equivalent loads were based on the expressions given in 
equation (6). The computational steps and results are detailed in Table 3. 

The relevant packet of input data cards is shown in Figure 3(c). 

Two computer runs were executed with the specified At-sets (b) and (c) . 
Table 6 displays their temperature results under the headings of Cases 12 and 13. 
Results are in excellent agreement, meaning that either At-set is satisfactory 
to render accurate solutions. The At-set (b) was then selected in studying the 
effect of the parameter P in the integration algorithm on temperature solutions. 

The integration algorithm available in the NTA for its transient thermal 
analysis is a modified Newmark Beta Method. It uses a fixed time-step size, 
but permits different time-step sizes to be changed at discrete times upon the 
user's specification. To avoid excessive matrix decompositions, the number of 
changes should be kept low. This algorithm is known to be an efficient one for 
large systems of linear equations, but it is uncertain for problems with strong 
nonlinearities. As accuracy and stability of a solution are closely related to 
the parameter (3 , a numerical investigation with different values of P was 
conducted. The values of p , including 0, 0.25, 1 and 0.55 (the default value), 
were executed. Except the case of using p = 0.5 (the central differences), 
whose temperature results are virtually identical to that of Case 13 (using p 
= 0.55) and, therefore, eliminated from the tabulation in Table 6. Case 16 
(using p = 1.0, the backward differences), yields comparable temperature results 
as that of Case 13. However, the solution instabilities caused both cases with 
P =0.25 and 0 (Euler integration) to terminate in their mid-courses. They 
are respectively labeled as Cases 17 and 18 and also included in Table 6. 

h = h(T r ) with T r = h(^l + T 200 ) 

When the average temperatures T r of the free stream temperature T 20 O anc * 
the surface temperatures of wall Ti are referenced by the temperature-dependent 
convective film coefficients h = h(Tr) , the MPC (Multipoint constraint) card is 
suited to provide such a relation. When T 100 = Tr i s assigned, the equivalent 
nonlinear thermal loads are functions of Ti00> meaning that Tioo(t) will be 
referenced by the nonlinear thermal load tables, if the modeling technique 
similar to that used in the two preceding cases is followed. Since Ti00 is 
constrained by an MPC, this modeling has violated one cardinal rule that must 
be strictly observed by a user in applying the nonlinear loads. The grid points 
to which the nonlinear loads are applied and the variables (either temperature 
or heat flux) on which they depend must not be variables eliminated by con- 
straints, i.e., they must be in the solution set. To circumvent this dilemma, 
a fictitious point GRID 150 was created for T 150 , and it was connected to T100 
by a super conductor, so that Tioo(t) = Ti5o(t). T 150 instead of T 100 was then 
referenced by. h(T r ), i.e., h(T r ) = h(Tioo) = MT 150 ), which, in turn, produced 
the equivalent nonlinear thermal loads. Figure 3(d) shows all relevant cards 
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and their relationships in detail, and the entered quantities in the nonlinear 
load tables are easily identifiable in comparison with those values given in 
Table 4. 

The two t-sets (b) and (c) were again selected to execute the computer 
runs. Results with the labeling of Cases 14 and 15 are shown in Table 7. 
Temperatures are still overall in good agreement with each other, but they no 
longer show the same quality as evidenced in previous cases. A scrutiny of the 
results reveals that relatively large temperature differences exist between 
these two temperature profiles in the range of time from 1.0 to 3.0 sec. This 
is a typical example showing the intricacies of this method , the interplays of 
the selected time-step sizes with the integration algorithm, and the effect on 
temperature results due to the presence of the MPC. 

Equation (4) shows that an inherent time lag error exists in applying the 
nonlinear thermal loads. They are computed for T of the (n+l)th time step 
based on temperature values at the nth and (n-l)th time steps. The larger the 
time-step size, the greater the deviation of the nonlinear thermal load from 
the desired value. When MPC is employed to relate a floating temperature (e.g., ; 
T]^) to a reference temperature (e.g., T 150 ) on which values of the nonlinear 
thermal loads depend indirectly, it introduces an additional step of time lag, 
and thus, an aggregated error for the nonlinear thermal load value. The errors 
are now induced by the time lag steps spanning from the (n+l)th to (n-2)th 
steps. Consequently, the selection of a finer time step size is appropriate 
in applying the nonlinear thermal loads, especially when an MPC is present 
to average two temperatures. 


OTHER APPLICATIONS AND RECOMMENDATIONS 


The modeling techniques of using the NLLEs to simulate the temperature- 
dependent convective film coefficients have been methodically illustrated. 

The demonstrated equivalent nonlinear thermal loads are all functions of 
temperature. If a simulation calls for the nonlinear thermal load to be a 
function of heat flux, a combination of the EPOINT and TF cards, in addition 
to the NOLINls and TABLEDls, have to be used. In all cases, values for the 
equivalent nonlinear thermal loads must be prepared by the user manually. This 
is an error-prone and elaborate process inherent in using the NLLEs. This 
deficiency, however, can be removed with only a relatively simple change made 
to the N0LIN1. It -will further broaden its capability greatly. The recommen- 
dation is to modify the computer code pertinent to S, the fifth field, in the 
N0LIN1 card so that it will accept any time- varying temperature, whether 
specified or computed values at some grid points in the model* The expanded 
capabilities will include, for instance, the same demonstrative problem with 
the hot-gas temperatures being a specified time-dependent temperature function. 
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CONCLUSIONS 


The nonlinear load elements have been illustrated to simulate nonlinear 
boundary conditions, in particular the temperature-dependent convective film 
coefficients in transient thermal analyses. The following remarks also serve 
as guidelines for those who attempt to use the NLLEs: 

(1) Prior to applying the nonlinear thermal loads to an NTA model in a 
transient thermal problem, a trial run on a simplified model is advisable. 

Only the linear thermal loads equal in magnitude to that of the nonlinear load 
set are necessary to apply to this model, and a relatively coarse time-step 
size may be used in the execution. The profile of the temperature results will 
provide a good indication on which the final selection of appropriate time-step 
sizes in runs with the real NLLEs depends. 

(2) Taking into consideration of the time lag errors, very small time-step 
sizes are always desirable. However, it is not necessary to select the same 
fine size uniformly over the entire range of time for integration. To maximize 
the processing efficiency, the selection of nonuniform time-step sizes according 
to the results of the trial run is a sensible and practical approach. Very fine 
time-step sizes should be reserved to concentrate at which the temperature 
profile has the largest slope in its rate change of temperature curve. 

(3) When MPC (also TF, if used) is present in the model for establishing a 
reference temperature that is referenced by the NLLEs, such as N0LIN1, a ficti- 
tious super thermal conductor has to be used to connect to that original grid 
point of referencing. A new grid point at the other end of the super thermal 
conductor is then to be referenced by any NLLEs so as not to violate the cardi- 
nal rule in applying the nonlinear loads. 

(4) For solution stability, which also affects accuracy, the proper selection 
of a value for the parameter (3 of the integration algorithm depends on the 
degree of nonlinearity in presence. However, only the range of 0 . 5 < |3 < 1 
should be considered. 
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TABLE 1 

DIMENSIONS AND THERMOPHYSICAL 
PROPERTIES OF THE COMPOSITE SLAB (PIPE) 


Material 

Stainless Steel (a) 

Insulation- Phenolic (b) 

Thermal Conductivity 
k (W/ cm-k) 

0.16258 

0.00502 

Thermal Capacitance 
PCp (J/cm^-k) 

4.029 

1.725 

Thickness 
Ax (cm) 

0.30 

0.40 

Cross-section 
A (cm^) 

1.0 

1.0 
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TABLE 2 


THE EQUIVALENT NONLINEAR LOADS FOR 
THE CASE OF h = 0.467 W/cm 2 -k 


T200 

N200 = 

— 

hAT200 

_ 

Ti 

-N]_ = hATi 

1422.0 

664. 

074 

293 

136.831 





589 

275.063 





867 

404.889 





1144 

534.248 





1422 

664.074 


TABLE 3 

THE EQUIVALENT NONLINEAR LOADS 
FOR THE CASE OF h = h(Ti) 


Tl 

h(Ti) 

-Nl(Ti) = hATx 

t 200 

N 200( T l) = hAT200 

293 

0.132 

38.676 

1422 

187.704 

589 

0.231 

136.059 



328.482 

867 

0.314 

272.238 



446.508 

1144 

0.393 

449.592 



558.846 

1422 

0.467 

664.074 



664.074 


131 





















TABLE 4 


THE EQUIVALENT NONLINEAR LOADS 
FOR THE CASE OF h = h(T ave ) 


t 200 

Tl 

T ave 

k( T ave-* 

N 200( T ave ) = hAT 200 

- N l< T ave> " hAT l 

14: 

12 

293 

857.5 

0.132 

187.704 

113.190 



589 

1005.5 

0.231 

328.482 

136.059 



867 

1144.5 

0.314 

446.508 

272.238 



1144 

1283.0. 

0.393 

.558.846 

449.592 



1422 

1422.0 

0.467 

664.074 

664.074 













TABLE 5 


A COMPARISON OF TEMPERATURE RESULTS 
FOR THE CASES WITH A CONSTANT h 



Case No. I 

Time 

(see) 

1 

2 

3 

10 

11 

0 

293.00 

293.00 

293.00 

293.00 

293.00 

0.25 

575.82 

580.41 

578.73 

582.08 

581.00 

0.50 

665.64 

668.29 

667.90 

669.82 

668.75 

0.75 

723.79 

725.67 

725.68 

727.06 

726.06 

1.00 

768.38 

769.88 

770.07 

771.19 

770.23 

1.25 

806.14 

807.42 

808.22 

809.19 

807.75 

1.50 

839.92 

841.08 

842.04 

842.93 

841.40 

1.75 

871.04 

872.11 

873.17 

873.99 

872.42 

2.00 

900.10 

901.10 

902.23 

903.00 

901.41 

3.00 

1000.39 

1001.21 

1003.56 

1004.18 

1002.48 

4.00 

1080. 22 

1080.90 

1083.37 

1083.86 

1082.23 

5.00 

1143.93 

1144.50 

1146.93 

1147.32 

1145.81 

6.00 

1194.83 

1195.29 

1197.61 

1197.92 

1196.53 

7.00 

1235. 52 

1235.90 

1238.05 

1238.29 

1237.04 

8.00 

1268.08 

1268.39 

1270.35 

1270.55 

1269.43 

HI 

(A) 

(C) 

(A) 

(B) 

- ... (c) _ 

CPU/I-0 

(min) 

0.125/1.170 

0.189/1.133 

0.154/1.309 

0.177/1.368 

0.198/1.256 
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TABLE 6 


A COMPARISON OF TEMPERATURE RESULTS 
FOR THE CASES WITH h = h^) 



Case No. 

Time 

(sec) 

12 

13 

16 

17 

18 

0 

293.00 

293.00 

293.00 





0.25 

401.14 

401.27 

401.06 





0.50 

445.89 

446.07 

445.82 





0.75 

478.54 

478.76 

478.61 





1.00 

505.58 

505.83 

505.82 





1.25 

529.69 

530.07 

530.17 





1.50 

552.17 

552.62 

552.80 





1.75 

573.65 

574.15 

574.40 





2.00 

594.44 

594.99 

595.28 





3.00 

673.30 

673.96 

674.38 





4.00 

746.00 

746.73 

747.09 . 





5.00 

813.37 

814.15 

814.46 





6.00 

875.67 

876.52 

876.80 





7.00 

929.48 

930.38 

930.66 





8.00 

977.03 

977.96 

.978.18 





At-sec 

(C) 

(B) 

(B) 

(B) 

(B) 

CPU/I-O 

(min) 

0.224/1.487 

0.200/1.606 

0. 192/1.615 

- 


- 

Parameter P 

0.55. 

0.55 

1.0 

0.25 

0 
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TABLE 7 


A COMPARISON OF TEMPERATURE RESULTS 
FOR THE CASES WITH h = h(T ave ) 


ec; 


Case No. 


14 


15 


0 

293.00 

293.00 

0.25 

350.94 

351.29 

0.50 

387.00 

387.63 

0. 75 

420.17 

420.99 

1.00 

453.04 

454.05 

1.25 

486.41 

488.16 

1.50 

521.83 

524.19 

1.75 

559.79 

562.77 

2.00 

600.52 

603.20 

3.00 

724.14 

725.32 

4.00 

832.04 

833.01 

5.00 

928.13 

928.48 

6.00 

1006.74 

1006.81 

7.00 

1075.18 

1075.07 


1135.05 

1134.79 

At-sec 

(B) 

(C) i 

CPU/I-0 

(min) 

0.192/1.661 

0.235/1.490 
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Figure 1 -vThe Composite Slab '(Pi-P.e.) and its. Finite . Element Model, 
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$ GRID 

POINTS 



GRID 

1 
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GRID 

2 


0.05 


GRID 

3 


0.10 


GRID 

A 


0.15 


GRID 

5 


0.20 


GRID 

6 


0.25 


GRID 

7 


0.30 


GRID 

8 


0.70 


GRID 
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-1.0 
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1 
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1 
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3 
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3 

A A 
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5 
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PROD 

101 

102 

1.0 
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7 
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0 

o 

• 
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S TIME 
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3 

5 
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+ TS2 
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Figure 2 — A Listing of the Input Data Deck of the Demonstrative Problem 
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S CUNVtOiIVE BOUNDARY IN CONTACT WITH HOT GAS 
CHbDY 1001 1002 PUINT 1 

+CY 200 

PHbUY 1002 1003 1.0 

i-1 AT 4 10 0 3 0.467 

(a) A constant h by the conventional modeling technique 


■ N 0 L INI 40 5 1 1 — i.o 1 

TABGED1 4051 

+TB1 293.0 136.831 589.0 275.063 867.0 

+TB2 1422.0 664.074 £NDT 

’NOG INI 405 1 1 664.074 1 

TA8GED1 4052 

+TB52 293.0 1.0 1450.0 1.0 ENDT 


1 4051 

+ 181 

404.889 1144.0 534.248 +T82 

1 40 52 

+ TB52 


(b) A constant h by the equivalent nonlinear thermal load approach 
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ENDT 
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= h(Tj_) by the equivalent nonlinear 

thermal load approach 
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(d) h = h(T r ) by the equivalent nonlinear thermal load approach 


Figure 3 - The Images of Packets of Data Cards for Simulating 
Various Convective Film Coefficients 
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NAS TRAN ANALYSIS OF HEAT- TRANSFER FLUID FILL PIPE FAILURES 


J. Ronald Winter 
Tennessee Eastman Company 



SUMMARY 

Failure analysis can be a very demanding endeavor that requires a great 
deal of interdisciplinary assistance if the correct solution to a particular 
problem is to be found. This paper presents an example that shows the 
difficulties one can encounter in such analyses and the advantage of the finite 
element method (NASTRAN) in assisting one in determining the true cause of a 
failure. In this example, cracks were developing along a pipe weld. After 
discarding several possible causes for the failures, it was finally determined 
that the problem was due to stress-corrosion-cracking associated with a rather 
unusual and novel environmental condition. 


INTRODUCTION 

In performing failure analysis an engineer often investigates numerous 
possible causes of a failure before he finally determines the real reason or 
reasons for the particular failure. This is usually due to either the lack of 
sufficient information or incorrect information. In many cases he must also 
either prove or disprove the reasons for the failure put forth by others. This 
situation is especially true in domestic industry and to a somewhat lesser 
extent in non-domestic establishments. In this paper, all of these important 
facts are presented rather than just the final solution. This should help new 
engineers in developing the logic for attacking such problems, and to be aware 
of the pitfalls and other oddities with which one must contend in order to 
obtain a successful solution. Hopefully, it will also make new engineers more 
aware of the need for assistance from other fields of engineering. Such inter- 
disciplinary assistance is often a necessity for solving certain problems. 


BACKGROUND 

At the time this particular problem was presented to the author, five 
previous failures had been encountered. The failures consisted of cracks in the 
fill pipe (stand pipe) of a closed (noncirculating) Dowtherm* (heat-transfer 


*Trademark of Dow Chemical Company. 
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fluid) system used to maintain the temperature of several spin blocks in a 
polyester spinning plant. The basic layout of the assembly which is relevant 
to this presentation is shown in Figure 1. The fill pipe and the location of 
the failure are shown in Figure 2. Detailed dimensions of the fill pipe are 
shown in Figure 3. The initial five failures consisted of cracks that mimicked 
bending fatigue failures. This had led several persons to believe that the 
failures were due to vibration induced by surrounding equipment or from some 
large vibrating dryers at another location in the building. These particular 
dryers had previously produced severe structural oscillations during transient 
conditions (start-up, shutdown, and coast down). This was the basic situation 
that existed when the author started the investigation. 

The need for an accurate solution to this problem increased dramatically 
as the number and frequency of failures increased due to the associated in- 
crease in lost production and the fire hazard associated with heat-transfer 
fluid leaks. Several minor fires were encountered before the problem was 
corrected. 


INITIAL ANALYSIS 

Although the failures did have the appearance of bending fatigue, I 
doubted this explanation because of the relatively large stiffness of the fill 
pipe assembly and the fact that any large oscillations that had been encounter- 
ed were in the Z direction (axial to the pipe) rather than in the X or Y 
(lateral) directions necessary to cause bending fatigue failures. But this 
remained to be proven. 

Since I was working on numerous other jobs at this time, the easiest and 
fastest way to perform the required vibration analysis was undertaken. This 
involved using a very simple NASTRAN finite element model* of the fill pipe. 
Using this method allowed me to easily account for the gusset stiffness and the 
tapered section of the flange. First, a normal mode (real eigenvalue) analysis 
was performed. This showed that the first natural frequency of this system was 
approximately 65 cps which, as expected, was far beyond the operating frequency 
of any of the vibrating equipment or electrical motors. The same model was 
then used to perform a forced response analysis using the transient oscilla- 
tions associated with the vibrating dryers as the forcing function. See 
Figure 4 for a typical response plot. This further indicated that vibration 
was not the problem. A frequency response analysis was also performed via 
NASTRAN which indicated the same results. With this data in hand, it was then 
concluded that the failures were not associated with vibration. This fact was 
forwarded to the production and maintenance personnel along with a request to 
save the next failure for careful inspection of the failure surfaces. All 
previous failures had been repaired by either cutting out the cracked region 
and rewelding or by installing a new section of pipe. 


*The model consisted of BAR elements for the pipe, plate elements for the 
gussets and a C0NM element for the heavy flange. See Figure 3. 
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Within a few weeks, two additional failures were encountered. One occur- 
red on a new assembly and one was a repeat failure. The repeat failure had 
been in service about three months. Most of the previous failures had been 
encountered in three to six month intervals. But there were several assemblies 
that had not failed even though they had been in operation for over a year . 
Thus, it became evident that a considerable amount of information was missing. 

As requested, the last two failures were held for me. They were taken to 
the materials laboratory where a micro-hardness test was performed. The micro- 
hardness test results indicated a substantial decrease in ductility in the weld 
region near the cracks. Due to an unusually heavy materials lab work load, no 
further metallurgical examinations were performed at this time. 

From the preparation of the initial NASTRAN model, I had become quite 
aware of the different materials of construction; i.e., stainless steel and 
carbon steel. The difference in the coefficients of thermal expansion for 
these materials would result in a differential thermal dilation (radial 
expansion) in the weld region. Such a situation results in the development of 
discontinuity stresses. Such discontinuity stress might at least partially 
explain the large hardness variations in the axial direction at the weld. Of 
course, some hardness variations were expected due to the residual stresses 
associated with welding. However, the variations were greater than anticipated 

To evaluate the discontinuity stresses in the material transition region, 

I decide^ £o apply the theory of beams on elastic foundations to this 
problem. ’ Two cases were employed in an attempt to bound the cause of the 
failures; i.e., one situation with conservative boundary conditions and the 
other with nonconservative boundary conditions. The first case consisted of 
assuming that the carbon section (flange) was infinitely stiff relative to the 
stainless steel pipe. These boundary conditions are shown below. The appli- 
cable equations are shown in Appendix A, Case A. 


CONSERVATIVE CASE 



Carbon Steel 


Is 

N 




X 2 


Analytical Constraints 
0, = 0 
Ar= Yi +V2 
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These relationships predicted stress levels far beyond the yield point 
of the material [>206.85 MPa (30,000 psi) ] . Since there was no evidence of 
such exceedingly high stress levels, this result was discarded. 

The second and hopefully nonconservative boundary conditions consisted of 
treating the weld region as the juncture of two infinitely long pieces of 
pipe, one made of stainless steel, the other of carbon steel. The boundary 
conditions for this situation are shown below. The applicable equations are 
presented in Appendix A, Case B. 


/ 

i 


i 


NON-CONSERVATIVE CASE 




Analytical Constraints 

P, =P 2 = P 
-M, = M 2 = M 
9, = e 2 

Ar = y, +y 2 


This analytical condition also predicted stress levels well beyond yield. 
From past experience, I knew that stress levels of this magnitude would pro- 
duce local deformation that one could easily feel and in most cases could be 
seen by the naked eye. However, a careful inspection of the samples did not 
reveal the evidence necessary to support these conditions. 

From this dilemma came the realization that the assumption of an abrupt 
change in material properties at the weld was actually incorrect. Proper 
treatment of the material properties along the fill pipe especially in the 
weld region should yield a more reasonable answer. This would, however, be a 
relatively time consuming, if not impossible, task to employ using the theory 
of beams on elastic foundations. But the finite element method would allow 
one to easily simulate materials variations and thus obtain reasonably good 
results even with a relatively coarse model. The more refined the model, the 
more accurate the solution. The axisymmetric nature of the stand pipe led to 
the development of two quarter models as shown in Figures 5 and 6 . 

The most important ingredients needed to obtain accurate analytical 
•results involved the determination of the material properties and a 
reasonable estimate of the wall thickness of the elements in the weld region. 
The thickness data was obtained by measuring several sliced sections from a 
typical pipe weld. An estimate of the material properties across the weld was 
obtained from some special reports supplied by our metallurgist. 

Both of the resulting models allowed one to vary the thickness and 
material properties as a function of the axial direction, Z or X, . (See 
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Appendix B for more model details.) In both models the following assumptions 
were made: 


1. The loads (thermal deformations) are axisymmetric with respect 
to the X. axis. 

2. The material properties vary along the weld region in the X. 
direction but are constant in the circumferential (hoop) 
direction. 

3 . Hookes law is obeyed . 

Model "a" consists of two straight sections of pipe, one 304 stainless 
steel, the other carbon steel. This situation, which is very similar to the 
second classical analysis (Case B) , resulted in thermally induced stress levels 
below the yield point of the materials. 

Model "b" is very similar to Model "a" but in this case the boundary 
conditions in the tapered section of the carbon steel were modified to simu- 
late the rigidity of the flange. This model very closely simulates the real 
boundary conditions. The stress levels were somewhat higher than those 
determined from Model "a," but still below yield. Close inspection of the 
weld region indicated that a much more refined model would be required to 
determine the peak stress levels at the beginning of the weld near the heat 
affected zone. A plot showing the thermally deformed weld region is presented 
in Figure 7 . 

To compensate for the relatively coarse grid network used in Models "a" 
and "b", a suitable stress concentration factor was determined. Application of 
this stress amplication factor indicated that the peak stress levels in the 
region of the failure would be between 137.9 MPa and 172.37 MPa (20,000 to 
25,000 psi) depending on the quality (roughness) of the weld. 

These results definitely showed that one must reasonably account for 
material property variation across a weld joining two different materials if 
there is to be any hope of obtaining a reasonable answer. However, these 
stress levels, which agreed with the deformation patterns determined from 
visual inspection of two failed pipes, left me with another dilemma. At these 
stress levels and only a few load cycles, failure should not occur. In addi- 
tion, the ductility should not have decreased substantially. This region 
would have had to be precracked during fabrication or there was something yet 
to be uncovered. 

At this point, I talked with, some of the maintenance personnel who were 
performing the field welding. They stated that occasionally when they tried to 
reweld the cracked region after proper grinding of the crack, the material would 
essentially crumble when they began welding. The material was seemingly very 
brittle, i.e., had lost its ductility. This was further evidence that some 
additional mechanism was contributing to the failures. 
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Fortunately, by this time the TEC Materials Laboratory was able to mount 
samples of a typical failure and obtain metallographic* photos. From these 
photos, the reason for the failures became obvious. It was stress-corrosion- 
cracking, 8003*4,5,6^ as evident in Figure 8. By this time failures of 
304 S.S. bourdon tubes and stainless steel rupture disks were also being 
encountered. Typical metallographic photos of these failures are shown in 
Figures 9 and 10. Thus, the puzzle began to fit together. We had determined 
that relatively high stress levels existed in the weld region. We were also 
aware that relatively large residual stresses could be present in the heat 
affected zone below the weld. Such stress levels in the presence of the proper 
chemical would readily explain the fill pipe failures. Usually, chlorides are 
the chief suspect for SCC in 300 series stainless steels, but no chloride 
source was immediately obvious. Thus we had to establish that chlorides were 
at fault and not some other chemical, then determine the source of the chlor- 
ides or other chemical. 

This led to a review of the procedure for manufacturing the spin blocks, 
the cleaning procedure and chemicals used for cleaning, and a check of the 
chloride concentration in the heat-transfer fluid being used in the spin 
block. 

This review was essentially fruitless. There was nothing that obviously 
stood out as a source of the chlorides. The chloride level in the heat- 
transfer fluid was found to be between 2 and 3 ppm. This concentration level 
was not nearly sufficient to cause the SCC being encountered. 

Again, we had a dilemma. At this time we contacted the Kodak Materials 
Laboratory at Kodak Park, Rochester, New York. Discussions with Kodak 
Material's engineers revealed that they had encountered corrosion problems 
in a closed system where too much moisture (H^O) was present. They had found 
that the moisture would move to the coolest part of the system where condensa- 
tion and revaporization would occur repeatedly. In our system such an internal 
reflux condition would serve to leach chlorides from the heat-transfer fluid 
and concentrate. them in the revaporization zone. 

Subsequent infra-red temperature measurements substantiated the reflux 
theory. As shown in Figure 11, the temperatures in the upper carbon steel 
portion of the fill pipes that had previously failed were quite cool in rela- 
tion to the heat-transfer fluid temperature. Thus, the failures in this 
relatively high temperature environment were attributed to stress-corrosion- 
cracking which was a result of : 

1. The welding of two materials with different coefficients of 
thermal expansion; 

2. The high operating temperature and the associated stress 
fields; and 

3. The presence of a mechanism to concentrate chlorides, i.e., a 
water/steam reflux cycle. 

*A standard polishing procedure with an acid etch, was employed. 
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The fill pipes that had not failed, did not show a temperature change in 
the weld to flange region. This indicated that these units did not contain 
sufficient moisture to establish the required water/steam reflux system for 
concentrating chlorides. 

From these results it was concluded that the problem could be corrected 
by either changing the flange material or by removing the moisture from the 
system. The latter corrective action was taken. 

Thus, through the cooperation and assistance from maintenance personnel 
(welders), production personnel, and material engineers at TEC and Kodak, I 
was finally able to establish the complete mechanism that was causing the 
failures. Such, interdisciplinary interactions are often a necessity in 
solving such unusual problems. 


CONCLUSION 

This analysis clearly shows the merits of the finite element method as 
available in NASTRAN. This capability allows the engineer to much more 
closely simulate real life situations and thus obtain more reasonable /accurate 
results. This is particularly helpful in failure analysis since conservative 
results do little more than cover up the real cause of a failure. Only in the 
design process should one consider conservative methods. But even in this 
case, one must be very cautious. 
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Figure 7 : NASTRAN Deformed Plot of Model b 



Figure 8: Photomicrograph Showing Stress Corrosion Craeking in the 304 

Stainless Steel Transition Region of the Fill Pipe 
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Figure 10: Additional Photomicrographs Showing Stress Corrosion Cracking 

in Other Fill Pipes and in the 304 Stainless Steel Pressure 
Relief Piping System 




Temperature (°F) 
Spin Block No. 



Figure 11: Thermal Profiles of Fill Pipes That Had Previously Failed 
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APPENDIX A 


Nomenclature : (See References 1 and 2.) 


P^^ = Shear load (lbf) 

LL = Moment (in- lbs) 

= Modulus of Elasticity (psi) 
u = Poisson's Ration 
r^ = Outside Radius of Pipe (inches) 


t^ = Wall Thickness (inches) 
T = Temperature (°F) 

0., = Slope 

y^ = Deflection (inches) 

I ± = bt 3 /12 = t 3 /12 


Case A: Conservative Case 

Stainless Steel Pipe to a Rigid Carbon Steel Block 

Assumption : Stainless steel pipe is infinitely long, thus semi-infinite 

beam on elastic foundation theory is applicable. 


i 

j 

I 

I 


CONSERVATIVE CASE 



Carbon Steel 


\ 

A 

A 






Analytical Constraints 
0 , =0 
Ar=y, +Vj 


Analytical Conditions and Data 
y 4 0 
©1 - 0 

T = 536°F = T 1 = T 2 


= 9.76 x 10 ^ in/in/°F 
~ 7.12 x 10 ^ in/in/°F 
r 1 = .73" 
t x = .137" 


Differential Radial Expansion 



Applicable Equations 


'ex 


'i =[ p i / 2 e i i i s i ] v - [y 2 e iv! 


e i ■ - [ V^iYi 2 ] V + [YWi] D 


ex 


Final Equations 

(1) = Ar 

which leads to 

j -•175P 1 + .716M-L = -14,400 

(2) B 1 = 0 

which leads to 
jP - 28^ = 8.18M 1 


Case B ; Non- Conservative Case 
Carbon Steel - Stainless Steel Pipe Junction 

Assumption: Both pipes are infinitely long, i.e., semi-infinite beams on 

an elastic foundation. 



NON-CONSERVATIVE CASE 




Analytical Constraints 

P, = P, = P 
-M, = M, = M 
0, =0, 

Ar = y, +y ? 
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Analytical Conditions and Data 


-M. 

= 

M = M 

a , 

— 

9.76 x 

10 

1 


2 

1 



-6 

P 1 

■ = 

P 2 = P 

°2 

= 

7.12 x 

10 

®1 

= 

<*2 

r l 

= 

.73" 


Ar 

= 

y l + y 2 

T 2 

= 

.783" 


T 

= 

536°F = T = T 2 

fc l 

= 

.135" 


AT, 

= 

at 2 

t 2 

= 

.240" 



Differential Radial Expansion 

Ar = r 1 d 1 (AT 1 ) - r^ATj = y ± + y 2 
Ar = ATCr^ - r^HT - T Q ) 


Applicable Equations 


y i = 
0 i = - 


V^iVi] d 8x - [v 2 Vi B i] 


Bx 

Bx 


J i=Y 


3(1 - u 2 )/r 2 t 2 


- -3* 


Bx 


e ' [jcos (Bx) - 

B gx = e Bx [sin (Bx)] 

C o =e' Bx 
Bx Q 

d r - e 
Bx 


cos ( Bx) - 
cos (Bx) ] 


Final Equations 

(1) Ar = y ± + y 2 

which leads to 
. 5424P - . 37 3M = 7200 

(2) 0 1 = 0 2 or © 2 - 0 2 = 0 

which leads to 
P = 33.04M 


in/in/°F 

in/in/°F 


sin (Bx) ] 
sin (Bx)] 
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A NASTRAN INVESTIGATION OF SIMULATED PROJECTILE DAMAGE EFFECTS 

ON A UH-lB TAIL BOOM MODEL 

Arnold T. Futterer 

U. S. Army Armament Research and Development Command 
Ballistic Research Labratory 

SUMMARY 

A NASTRAN model of a UH-lB tail boom that had been designed for a another 
project was used to investigate the effect on structural integrity of 
simulated projectile damage. Elements representing skin, and sections of 
stringers, longerons and bulkheads were systematically deleted to represent 
projectile damage. The structure was loaded in a manner to represent the 

flight loads that would be imposed on the tail boom at a 130 knot cruise. The 
deflection of four points on the rear of the tail boom relative to the 
position of these points for the unloaded, undamaged condition of the tail 
boom was used as a measure of the loss of structural rigidity. The same 
procedure was then used with the material properties of the aluminum alloys 
replaced with the material properties of T300/5208 high strength 
graphite/epoxy fibrous composite material, (0,±45,90) g for the skin and 
(0,±45) g for the longerons, stringers, and bulk heads. 

INTRODUCTION 

This investigation had a two-fold objective : 

1. To determine the effect on the structural integrity of the UH-lB tail 
boom caused by threat projectile damage. 

2. To estimate the effect of composite materials on the stiffness of the 
tail boom. 

The model of the UH-lB tail boom used in the analysis was originally 

prepared under contract by Kamen AviDyne^ (KAD) for the Ballistic Research 
Laboratory. The model consisted of beams representing sections of the 
stringers, bulk heads, and longerons and thin plates representing the skin. 
The KAD report describes the model in good detail. Figures 1, 2, and 3 taken 

from the KAD report illustrate the NASTRAN^ model developed by KAD and give 
the numbering schemes for the grid points, beam elements and plate elements, 
respectively. The skin is made of 2024 T3 aluminum alloy with a modulus of 

elasticity of 7.31x10^ MPa (10.6 x 10^ psi) and a mass density equal to 271.15 

kg sec^/m^ (0.00025 lb sec^/in^). The stringers, bulk heads and longerons are 
made of 7075 T6 aluminum alloy with a modulus of elasticity of 7.10xl0 2 * 4 MPa 

(10.3 x 10^ psi) and a mass density of 271.15 kg sec^/m^ (.00025 lb 

sec^/in ). For further detail on the assumptions that went into preparing the 


161 







FIGURE 3. UH-1B TAIL BOOM MODEL PLATE ELEMENT IDENTIFICATION 



model it is recommended that a copy of the report be obtained from Defense 
Technical Information Center. 


PROCEDURE 

The investigation was accomplished by using the NASTRAN code, a complex 
finite element method computer program. In all 35 NASTRAN calculations were 
made, 17 using the aluminum alloy construction and 18 using the same structure 
but using the material properties of T300/5208 high strength, graphite-epoxy 
fibrous composite. This paper describes the results from 18 calculations (9 
aluminum alloy construction and 9 composite construction) The relative 
displacements, compared to tne non-damaged, non-loaded structure, of the four 

points at the base of the tail fin with the front of the tail boom anchored 

securely were used as a measure of the deterioration of the structural 
integrity of the tail boom due to projectile damage. The flight loads that 
would be imposed on a helicopter cruising at 130 knots were simulated by 
loading the structure to simulate the loads and torques that the rotor thrust 
and elevator loading would cause at cruise velocity. The assumption was made 

that a large hole or tear in a skin panel or a break in a longeron, stringer 

or bulkhead section would destroy the structural integrity of the element 
representing that skin panel or section. In the model that damaged element 
was then deleted. The investigation was conducted by systematically deleting 
plate and beam elements to simulate greater damage. Table I gives the 

nomenclature of the damage configurations that were investigated. Damage to 

TABLE I. NOMENCLATURE FOR DAMAGE CONFIGURATIONS 
Nomenclature Elements deleted 


0 No elements deleted 

1 419 


2 

419, 

141, 

418 





3 

431, 

153, 

419 





4 

431, 

153, 

419, 

141, 

418 



5 

419, 

154, 

431, 

165, 

430, 

153, 

418, 141 

6 

419, 

154, 

431, 

165, 

430, 

153, 

0 

418, 141, 164, 140 

7 

458 







8 

458, 

222 
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the left and to the right side of the tail boom were studied because the right 
side has thicker skin than the left. Furthermore, the longitudinal strains 
were generally compressive on the right side and tensile on the left side. 
This paper discusses the results of calculations for damage to the right side 
of the helicopter. The numbers of the elements deleted refer to Figures 2 and 
3. As may be noted, the 100 and 200 series numbers refer to the beam elements 
and the 400 series to the rectangular plate elements. 


Tables II, and III are similar in construction. They give the 
displacements of the points at the rear of the tail boom of the loaded, 
undamaged tail boom and the loaded tail boom with simulated damage relative to 
the undamaged, unloaded state of the tail boom. Nomenclature refers to Table 
I where the damage configurations are set forth. Material lost refers to the 
mass of the deleted elements. "Direction" ,X, Y, and Z, gives the displacement 
of the grid points in the three coordinate directions given on Figure 1 and 
"R", which is the square root of the sum of the squares of the three 
coordinate displacements, gives the total displacement of the grid points 
specified. The displacements are given in inches and millimeters. The 
displacement values are given in exponential format, i.e., 2.05-2 means 2.05 x 


—2 

10 . Table II contains the results of the calculations with the tail boom 
constructed of aluminum alloy and the damage is to the forward right side of 
the helicopter. Figure 4 is a graph of deflection of points on the rear of 
the tail boom versus mass of aluminum alloy removed from the forward right 
side of the tail boom. As may be noted, the deflection is quite linear with 
mass removed until about 1 kg and then further removal causes the 
displacements to become non-linear suggesting a more rapid approach to failure 
with further loss of material. 


After calculating the displacements for the various damage configurations 
with the tail boom constructed of aluminum alloys, the calculations were 
repeated using the material properties of T300/5208 which is a high strength, 
graphite-epoxy fibrous composite. The skin plates were assumed to be 
constructed of [0,±45,90] g layered , composite and the beam elements of [0,±45] g 

layered composite. T300/5208 was recommended as being high strength and 
considerably less expensive than the ultrarhigh modulus graphite-epoxy. Since 
the composites have less strength in compression than in tension, the material 
moduli of elasticity for both compression and tension were used in the 

calculations. For damage on the left side the tensile moduli, 5.59 x 10^ MPa 

(8.11 x 10^ psi) for the plate elements and 6.50 x 10^ MPa (9.42 x 10^ psi) 
for the beam elements, were used. For damage on the right side the 

compressive moduli, 5. 38 x 10^ MPa (7.81 x 10^ psi) for the plate elements and 

6.43 x 10^ MPa (9.32 x 10^ psi) for beam elements, were used. A mass density 

of 162.69 kg sec^/m^ (.00015 lb sec^/in.^) was used for the beam and plate 
elements. This paper reports on the more extreme of the two cases, the right 
side damage and using the lower compressive moduli in the calculations. The 
total displacements for the various damage configurations for damage done to 
the forward right side of a tail boom constructed of T300/5208 composite are 
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TABLE II. DAMAGE TO THE RIGHT SIDE 


Nomen- 

Material 

Direc- 



clature 

Lost 

tion 


67 


lb/kR 


in. 

mm 

0 

.0 

X 

2.18-2 

5.54-1 


.0 

Y 

3.30-1 

8.38+0 



Z 

-2.22-1 

-5.64+0 



R 

3.98-1 

1.01+1 

1 

.71 

X 

2.17-2 

5.51-1 


.32 

Y 

3.36-1 

8.53+0 



Z 

-2.24-1 

-5.69+0 



R 

4.04-1 

1.03+1 

2 

1.46 

X 

2.14-2 

5.44-1 


. 66 

Y 

3.60-1 

9.14+0 



Z 

-2.36-1 

-5.99+0 



R 

4.30-1 

1.09+1 

3 

1.46 

X 

2.17-2 

5.51-1 


. 66 

Y 

3.43-1 

8.71+0 



Z 

-2.28-1 

-5.79+0 



R 

4.12-1 

1.05+1 

4 

2.33 

X 

2.13-2 

5.41-1 


1.06 

Y 

3.64-1 

9.25+0 



Z 

-2.40-1 

-6.40+0 



R 

4.36-1 

1.11+1 

5 

3.00 

X 

2.12-2 

5.38-1 


1.36 

Y 

3.78-1 

9.60+0 



Z 

-2.51-1 

-6.38+0 



R 

4.54-1 

1.15+1 


’ HELICOPTER TAIL BOOM CONSTRUCTED OF ALUMINUM ALLOY. 
Grid Point Displacements 

70 71 72 

in. mm in. mm in. mm 


2.45-3 

6.22-2 

7.52-3 

3.30-1 

8.38+0 

3.23-1 

-2.47-1 

-6.27+0 

-2.46-1 

4.12-1 

1.05+1 

4.06-1 

2.02-3 

5.13-2 

7.35-3 

3.36-1 

8.53+0 

3.29-1 

-2.50-1 

-6.35+0 

-2.50-1 

4.19-1 

1.06+1 

4.13-1 

2.21-3 

5.61-2 

6.49-3 

3.60-1 

9.14+0 

3.53-1 

-2.64-1 

-6.71+0 

-2.62-1 

4.46-1 

1.13+1 

4.40-1 

1.55-3 

3.89-2 

7.29-3 

3.43-1 

8.71+0 

3.34-1 

-2.55-1 

-6.48+0 

-2.53-1 

4.27-1 

1.08+1 

4.19-1 

-1.99-3 

-5.05-2 

6.63-3 

3.64-1 

9.25+0 

3.54-1 

-2.70-1 

-6.86+0 

-2.66-1 

4.53-1 

1.15+1 

4.43-1 

-1.24-3 

-3.15-2 

6.04-3 

3.78-1 

9.60+0 

3.69-1 

-2.82-1 

-7.16+0 

-2.77-1 

4.72-1 

1.20+1 

4.61-1 


1.91-1 

-1.23-2 

-3.12-1 

8.204-0 

3.23-1 

8.20+0 

-6.25+0 

-2.70-1 

-6.86+0 

1.03+1 

4.21-1 

1.07+1 

1.87-1 

-1.27-2 

-3.23-1 

8.36+0 

3.29-1 

8.36+0 

-6.35+0 

-2.74-1 

-6.96+0 

1.05+1 

4.28-1 

1.08+1 

1.65-1 

-1.49-2 

-3.78-1 

8.97+0 

3.53-1 

8.97+0 

-6.65+0 

-2.88-1 

-7.32+0 

1.12+1 

4.56-1 

1.16+1 

1.85-1 

-1.31-2 

-3.33-1 

8.48+0 

3.34-1 

8.48+0 

-6.43+0 

-2.79-1 

-7.09+0 

1 . 06+1 

4.35-1 

1.10+1 

1.68--1 

-1.50-2 

-3.81-1 

8.99+0 

3.54-1 

8.99+0 

-6.76+0 

-2.94-1 

-7.47+0 

1.13+1 

4.60-1 

1.17+1 

1.53-1 

-1.64-2 

-4.11-1 

9.37+0 

3.69-1 

9.37+0 

-7.04+0 

-3.06-1 

t7.77+0 

1.17+1 

4.79-1 

1.22+1 



TABLE II CONT. DAMAGE TO THE RIGHT SIDE OF HELICOPTER TAIL BOOM CONSTRUCTED OF ALUMINUM ALLOY 


Nomen- 

clature 

Material 

Lost 

lb/kg 

Direc- 

tion 


67 

Grid Point Displacements 

70 71 



72 

in. 

mm 

in. 

mm 

in. 

mm 

in. 

mm 

6 

3.90 

X 

2.12-2 

5.38-1 

-3.25-3 

-8.25-2 

4.45-3 

1.13-1 

-1.96-2 

-4.98-1 


1.77 

Y 

4.08-1 

1.04+1 

4.08-1 

1.04+1 

4.02-1 

1.02+1 

4.02-1 

1.02+1 



Z 

-2.73-1 

-6.93+0 

-3.04-1 

-7.72+0 

-3.01-1 

-7.65+0 

-3.31-1 

-8.41+0 



R 

4.91-1 

1.25+1 

5.09-1 

1.29+1 

5.02-1 

1.28+1 

5.21-1 

1.32+1 

7 

3.18 

X 

2.62-2 

6.65-1 

-6.38-4 

-1.62-2 

8.04-3 

2.04-1 

-1.07-2 

2.72-1 


1.44 

Y 

3.51-1 

8.92+0 

3.51-1 

8.92+0 

3.03-1 

7.70+0 

3.03-1 

7.70+0 



Z 

-2.47-1 

-6.27+0 

-3.33-1 

-8.46+0 

-2.74-1 

-6.96+0 

-3.48-1 

-8.84+0 



R 

4.30-1 

1.09+1 

4.84-1 

1.23+1 

4.09-1 

1 . 04+1 

4.62-1 

1.17+1 

8 

4.25 

X 

2.24-2 

5.69-1 

4.42-4 

1.12-2 

7.77-3 

1.97-1 

-1.10-2 

-2.79-1 


1.93 

Y 

3.53-1 

8.97+0 

3.53-1 

8.97+0 

3.03-1 

7.70+0 

3.03-1 

7.70+0 



Z 

-2.48-1 

-6.30+0 

-3.37-1 

-8.56+0 

-2.76-1 

-7.01+0 

-3.52-1 

-8.94+0 



R 

4.32-1 

1.10+1 

4.88-1 

1.24+1 

4.10-1 

1.04+1 

4.65-1 

1.18+1 




TOTAL DISPLACEMENT OF GRID POINTS 


0 1 2 3 libs) 4 



DAMAGE LOSS OF STRUCTURAL WEIGHT (kg) '7 

i 

FIGURE 4. GRID POINT DISPLACEMENT VS. DAMAGE LOSS OF STRUCTURAL WEIGHT-ALUMINUM ALLOY CONSTRUCTION 



compiled in Table III. Figure 5 is a graph of the relative displacement of 
points on the rear of the tail boom versus composite material lost through 
simulated projectile damage to the forward right side. As is to be expected, 
the non-linear behavior of the curve shows up at lower mass loss due to the 
lower strength of the material for a given cross-section. 

A further consideration is that the tail boom constructed of aluminum 
alloy had a structural weight of 57.7 kg (127.1 lb) and in the model 
additional non-structural weight amounting to 24.0 kg (52.9 lb) was 
distributed along the length of the tail boom. Replacement by the lighter 
weight composite would result in a weight saving of from 23.1 kg (50.9 lb) up 
to 32.7 kg (72.0 lb) depending to what extent the aluminum alloy could be 
replaced by the composite material. This reduction in weight could make the 
helicopter more maneuverable or able to carry a larger payload. A further 
consideration is whether use of the more expensive ultra-modulus composite 
with a modulus almost twice that of T300/5208 and about 50% greater than the 
aluminum alloy might be warranted. The density of the ultra-modulus 
graphite/epoxy is only about two percent greater than that of T300/5208. 


CONCLUDING REMARKS 

The study completed here leaves many questions while providing some 
answers. Subjects of further study should be to what extent: the loss of 
structural rigidity of the tail boom can be tolerated; does a composite 
react to projectile damage better or worse than the aluminum alloy. Can the 
price differential of the ultra-high modulus composite be tolerated in the 
construction of the tail boom taking . into consideration its much better 
material qualities? 
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TABLE III. DAMAGE TO THE RIGHT SIDE OF HELICOPTER TAIL BOOM CONSTRUCTED OF T300/5208 COMPOSITE. 


Nomen- Material Direc- 
clature Lost tion 
lb/kg 

0 .OX 

.0 Y 


1 .43 

.19 



in. 

mm 

in. 

mm 

in. 

mm 

in. 

mm 

2.78-2 

7.06-1 

2.83-3 

7.06-2 

9.36-3 

2.45-1 

-1.56-2 

-3.96-1 

4.22-1 

1.07+1 

4.22-1 

1.07+1 

4.12-1 

1.05+1 

4.12-1 

1.05+1 

-2.79-1 

-7.09+0 

-3.17-1 

-8.05+0 

-3.11-1 

-7.90+0 

-3.46-1 

-8.79+0 

5.06-1 

1.29+1 

5.28-1 

1.34+1 

5.16-1 

1.31+1 

5.38-1 

1.37+1 

2.77-2 

7.04-1 

2.28-3 

5.79-2 

9.17-3 

2.33-1 

-1.61-2 

-4.09-1 

4.29-1 

1.09+1 

4.29-1 

1.09+1 

4.19-1 

1.06+1 

4.19-1 

1.06+1 

-2.82-1 

-7.16+0 

-3.21-1 

-8.15+0 

-3.15-1 

-8.00+0 

-3.50-1 

-8.89+0 

5.13-1 

1.30+1 

5.35-1 

1.36+1 

5.24-1 

1.33+1 

5.46-1 

1.39+1 

2.73-2 

6.93-1 

8.68-5 

2.20-3 

8.14-3 

2.07-1 

-1.87-2 

-4.75-1 

4.59-1 

1.17+1 

4.59-1 

1.17+1 

4.49-1 

1.14+1 

4.49-1 

1.14+1 

-2.97-1 

-7.54+0 

-3.37-1 

-8.56+0 

-3.30-1 

-8.38+0 

-3.68-1 

-9.35+0 

5.47-1 

1.39+1 

5.69-1 

1.45+1 

5.57-1 

1.41+1 

5.81-1 

1.48+1 

2.76-2 

7.01-1 

1.72-3 

4.37-2 

9.12-3 

2.32-1 

-1.66-2 

-4.22-1 

4.37-1 

1.11+1 

4.37-1 

1.11+1 

4.25-1 

1 . 08+1 

4.25-1 

1.08+1 

-2.87-1 

-7.29+0 

-3.28-1 

-8.33+0 

-3.19-1 

-8. 10+0 

-3.57-1 

-9.07+0 

5.23-1 

1.33+1 

5.46-1 

1.39+1 

5.31-1 

1 . 35+1 

5.55-1 

1.41+1 

2.72-2 

6.91-1 

-3.99-4 

-1.01-2 

8.37-3 

2.13-1 

-1.88-2 

-4.78-1 

4.63-1 

1.18+1 

4.63-1 

1.18+1 

4.50-1 

1.14+1 

4.50-1 

1.14+1 

-3.02-1 

-7.67+0 

-3.46-1 

-8.79+0 

-3.34-1 

-8.48+0 

-3.75-1 

-9.53+0 

5.53-1 

1.40+1 

5.78-1 

1.47+1 

5.60-1 

1.42+1 

5.86-1 

1.49+1 

2.71-2 

6.88-1 

-1.67-3 

-4.24-2 

7.68-3 

1.98-1 

-2.05-2 

-5.21-1 

4.80-1 

1.22+1 

4.80-1 

1.22+1 

4.67-1 

1.19+1 

4.67-1 

1.19+1 

-3.14-1 

-7.98+0 

-3.60-1 

-9.14+0 

-3.48-1 

-8.84+0 

-3.97-1 

-1.01+1 

5.74-1 

1.46+1 

6.00-1 

1.52+1 

5.82-1 

1.48+1 

6.13-1 

1.56+1 




Table III CONT. DAMAGE TO THE RIGHT SIDE OF HELICOPTER TAIL BOOM CONSTRUCTED OF T300/5208 OOMPOSITE 


Nomen- 

clature 

Material 

Lost 

lb/kg 

Direc- 

tion 



Grid Point Displacements 





67 

70 

71 



72 

in. 

mm 

in. 

mm 

in. 

mm 

in. 

mm 

6 

2.34 

X 

2.70-2 

6.86-1 

-4.58-3 

-1.16-1 

5.39-3 

1.51-1 

-2.52-2 

-6.40-1 


1.06 

Y 

5.24-1 

1.33+1 

5.24-1 

1.33+1 

5.15-1 

1.31+1 

5.15-1 

1.31+1 



Z 

-3.46-1 

-8.79+0 

-3.93-1 

-9.98+0 

-3.83-1 

-9.73+0 

-4.26-1 

-1.08+1 



R 

6.28-1 

1.60+1 

6.55-1 

1.61+1 

6.42-1 

1.63+1 

6.68-1 

1.70+1 

7 

1.91 

X 

2.87-2 

7.29-1 

-3.00-5 

-7.62-4 

1.00-2 

2.54-1 

-1.39-2 

-3.53-1 


.87 

Y 

4.46-1 

1.13+1 

4.46-1 

1.13+1 

3.89-1 

9.88+0 

3.89-1 

9.88+0 



Z 

-3.07-1 

-7.80+0 

-4.13-1 

-1.05+1 

-3.42-1 

-8.69+0 

-4.32-1 

-1.10+1 



R 

5.41-1 

1.37+1 

6.08-1 

1.54+1 

5.18-1 

1.32+1 

5.81-1 

1.48+1 

8 

2.55 

X 

2.85-2 

7.24-1 

-5.14-4 

-1.32-2 

9.64-3 

2.45-1 

-1.44-2 

-3.66-1 


1.16 

Y 

4.48-1 

1.14+1 

4.48-1 

1.14+1 

3.90-1 

9.91+0 

3.90-1 

9.91+0 



Z 

-3.10-1 

-7.87+0 

-4.17-1 

-1.06+1 

-3.45-1 

-8.76+0 

-4.37-1 

-1.11+1 



R 

5.45-1 

1.38+1 

6.12-1 

1.55+1 

5.21-1 

1.32+1 

5.85-1 

1.49+1 





FINITE CIRCULAR PLATE ON ELASTIC FOUNDATION CENTRALLY LOADED 
BY RIGID SPHERICAL INDENTER 
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This paper discusses the analytical solution of a finite circular plate 
on an elastic foundation centrally loaded by the rigid indenter; and the 
procedure to use NASTRAN as a subroutine to iteratively converge to this 
solution numerically. 


INTRODUCTION 


In contact recording applications, where a thin flexible spinning disk 
backed by an elastic pad is penetrated by a read/write head, it is 
important to understand the head penetration, the contact area, and the 
pressure distribution. A computer program was written to iteratively 
converge to the solution using MSC/NASTRAN version 60* (hitherto called 
NASTRAN) as a subroutine. Results were erratic. It was found that 
extremely fine grid was needed around the contact area to converge to 
the right solution. 

In order to calculate the approximate contact area, the model was 
simplified to a finite circular plate on an elastic foundation centrally 
loaded by the head. The approximation is fairly valid because the 
flexural rigidity of the plate is extremely low and the reaction forces 
at the center of the spinning disk is about 2% of the total load. The 
effect of most of the load, thus, is to cause local deformation of the 
plate between the head and the pad. An analytical solution for this 
case was obtained. The iterative procedure to converge to this solution 
by using NASTRAN as a subroutine to verify the analytical solution was 
also defined. 


MacNeal-Schwendler Corp version of NAsa STRuctural ANalysis. 
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ANALYTICAL SOLUTION 


Consider in figure 1 the spherical indenter with radius being pressed 

into the circular plate with center at the origin, the plate resting on 
an elastic foundation with elastic spring constant k. It is of interest 
to determine the inside contact radius OA and the pressure distribution 
on the indenter over that contact. Notice, as shown in figure 1, that 
the outside contact radius is OB and the plate physically separates from 
the elastic foundation for radii larger than OB. Behavior of this 
"free" plate can also be modeled easily by matching boundary conditions 
at B, but it was not necessary in our problem. Instead we were seeking 
the solution to the problem "what is the outside radius R q of the plate 

such that vertical displacement W=0 at R^?" 

The behavior of the plate on elastic foundation can be described (ref. 1) 
by equation 


2 2 t r 2 2 

±L +11 +1 iL 1 1 l_w l aw l a_w 

. 2 + r 8r + 2 „.2 . 2 r 3r 2 ao 2 

L 3r r 30 J L 3r r 30 

For our case, in which the indenter presses on the center of the plate, 
there is a radial symmetry, and the solution is independent of angle 0, 
we can rewrite the above equation as: 


q-kW 

D 



+ 


1 d 
r dr 


d 2 W 1 dW 

72 r 37 

dr 


q-kW 

D 


( 1 ) 


Where W = vertical displacement 

q = distributed lateral load 
D = flexural rigidity of the plate. 


Eh 

12(1-7) 


h = thickness of the plate 
E = Young’s modules of the plate 
v = Poisson' s ratio of the plate 

k = spring rate/unit area of the elastic foundation. 
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In seeking the solution for above, we divide the plate into two sections, 
as shown in figure 2. The inside (section 2) portion conforms to the 
sphere up to contact radius R^. The outside portion of the plate must 

also satisfy differential equation (1) with appropriate matched boundary 
conditions . 


OUTSIDE PORTION OF THE DISK 


Consider now the outside portion of the disk (section 1 in figure 2) . 

If q = constant, such as gravity load, the particular solution of 
equation (1) is simply Therefore, we now seek the complementary 
solution (in absence of q) of 


2 2 2 
cT Id dTW 1 d W 

,2 r dr ,2 r dr 

dr J dr 


Let us define Z 



-kW 

D 


5 


r_ 

z 


W R 

n = s 


( 2 ) 


where R is the radius of the indenter. 
s 

Equation 2 can then be rewritten as 


r£ + I 

d 

' d 2 n , 1 

an 

.dS 2 5 

d? 

. 8? 2 ? 

95 


This is a linear differential equation of the 4th order. The solution 
can be represented in terms of Bessel's function. However, it is 
presented below in terms of series solution for ease of computation. 
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The general solution is given by 


n = A 1 «1 + A 2 «2 + A 3 <$ 3 + A 4 6 4 


( 3 ) 


where A. . . A, are constant, 
1 4 


IW1. 


<5, = 1- 


(4) 


2 2 2 2 2 2 
274 27476787 


2 _6 plO 

<5„ = r - 


C5) 


2 2 2 2 2 2 
476^ 47678710 


Coefficients in above two series can be generated by recursive relationship 


U , 

n _ n~4 

n y 2 % , ox 2 


(n ) (n-2) 


6 3 = 6 1 log £ + & 3 

64 = 6 1 log C + A 4 jJ 

where 

A = _J_ 25 C 8 ' (6) 

3 128 “ 1769472 

A = _J_ C 6 1540 C 10 (7) 

4 3456 " 442368 
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and coefficients in 6 and 7 can be recursively generated by 


, _ -1 . , 4. (n) (n-1) (n-2) (-1) 

n 2 , ' .2 n-4 0 2.2,2 2 

n (n-2) L 2.4.6. ...n 


n/ 4 


Four constant A^..A^ can be thus evaluated by following boundary conditions. 

dW r 


1 At r = R. 


Cl 


dr 


... conforms to the suhere 


R 


Therefore, at E, = 


Cl 


dn 

d£ 


dW 

d£ 


dW 

dr 


dr 

d? 


s_ . £ 
2 


fl.gi 

R 


- ? i 


At r = R, 


Cl 


d 2 W 

dr 2 


1 

R 


. . conforms to sphere 


d 2 n _ 

d € 2 


= -1 
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3 


Moment 


\ 


iortion of the disk can be obtained 


(. 



SOLUTION OF INSIDE PORTION OF DISK 


The solution for the inside portion of the disk (section 2, figure 2) 
that conforms to the sphere satisfies 


d 2 W 


, 2 R 
dr s 


and, therefore. 


w - »! + 2iT < R c? - r2) 

S 


( 8 ) 


where VL = deflection of the disk at r = R^ obtained from equation (3) 


PRESSURE DISTRIBUTION 


Total force F is given by 


F = 


m 

J 


R 

. o 


(kW) 2irr . dr 


where W = W i + 2R (R CI ° <r<R CI 

s 


A 

R 


A 1 S 1 + A 2 5 2 + A 3 5 3 + A 4 J 4 


E CI <r<K o 


Pressure distribution under the head is given by equation (8) multiplied 
by k. However, at radius R CI there is additional shear force required to 

keep the outside portion of the disk in equilibrium, and is given by 


Q = D 




1_ dW 

2 dr 
r 
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NORMALIZING 


Solution to the equations are functions of R g , F, £, R q . However, if 

the disk is very weak and outside radius R q is taken to be that value 

for which W=0, then from equations (3 through 8), it can be shown that 

FR , R A , OR £ , qR £^ are all normalized functions of 

s * o' x s ’ ■ s 

R ci/r 

Figures 3 and 4 show some of the relationships. 

NASTRAN MODEL DESCRIPTION 


The application of the static analysis of MSC/NASTRAN to the problem in 
figure 1 is demonstrated in this section. 

The elastic foundation in figure 1 is first discretized into the "gap 
scalar spring" of finite length corresponding to grid points on the 
plate. The "gap scalar spring" is defined as a linear spring that can 
be compressed only. An extension of this spring will produce no spring 
force. This spring force will then be used to generate the FORCE card 
for the NASTRAN program. Secondly, an isoparametric bending element 
with transverse shear is used to model the plate between the indenter 
and the elastic foundation. As shown in figure 5, to save computer 
time, only a section of this circular plate was modeled in this work. A 
symmetric condition is used along the edge boundaries. Since the spherical 
indenter can be treated as a rigid body in this problem, its contour is 
computed and stored in [C] . This [C] is used as an enforced displacement 
to the plate model within the contact area. 

The radius R of the contact area between the indenter and the plate is 

L* X 

essentially the solution we need. In the following part, an iterative 
technique was developed to compute the R^, contact radius, and the load 

distribution over the plate. 

Final solution to the problem must meet the following conditions: 

1 There should be no geometric interference between the deformed 

plate and the contour of the indenter. The interference example is 
as shown in figure 6, where the plate deformation pattern IW]^ 

interferes with the indenter contour [C] . The smallest R^ without 

a geometric interference is the solution to this problem. 
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2 The s umma tion of the distributed load, or the "gap scalar spring" 
force, over the plate be equal to the total force F. 

£ [q] [area] = F (9) 

In order to satisfy the above two conditions, an iterative type of 
algorithm was developed to solve the problem by utilizing the CALL 
NASTRAN technique. In this case, the whole NASTRAN program becomes a 
subroutine that can be implemented by any user-developed program. The 
NASTRAN program can thus be invoked by using the CALL statement in the 
user's program just like using any other conventional subroutines. 

To start this iteration process, an initial guess of R CI is required as 

an input to the program. With this assumed R , the static contact area 

between the indenter and the plate can be defined and set equal to the 
indenter contour by using the enforced displacement card, SPC. Initial 
plate deflections are obtained by using NASTRAN. With these deflections 
and the elastic foundation stiffness k, forces on the grid can be calculated 
satisfying condition (2) above. This [q], as the FORCE cards, along 
with the previously defined SPC cards, are then put into the NASTRAN 
program to compute the plate displacement [W] . This plate displacement 
is not supposed to interfere with the indenter contour IC] . If this 
geometric constraint is not satisfied, a larger R CI will be assumed next 

and repeat the above steps until a R^ can be found that satisfies this 

geometric constraint. This final R^ is the solution we are trying to 

obtain. Thus, the static contact area and the final distributed load 
can be computed accordingly without any difficulty. 

Figure 7 shows the simplified flow diagram of this iteration process. 

Some additional explanations and programming considerations are made 
below. 

1 Since the NASTRAN program is being called and invoked many times, 
NASTRAN CHECKPOINT and RESTART features have been used. 

2 The CALL NASTRAN statement actually includes the following functions: 

a Generate the NASTRAN executive control deck, including CHECKPOINT 
dictionary for RESTART use. Case Control Decks and the Bulk 
Data Card. 

b Invoke the NASTRAN program by using a CALL statement. The 
NASTRAN program is invoked throughout the entry point of its 
load module. This entry point name is "NASTRAN" in our case. 
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c Recover the displacement vector computed by NASTRAN as the 
output from this GALL NASTRAN operation. The module OUTPUT4 
was used to modify the DMAP sequence to get the desired output 
from NASTRAN. 


RESULTS 


The iteration scheme defined in the NASTRAN model gave results quite 

close to those obtained analytically. The graphs are shown in figures 3 

and 4. Given force F and the radius R of the indenter, contact radius 

s 

and the pressure distribution can be obtained from figures 3 and 4 
respectively. 
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SPHERICAL INDENTER WITH PLATE ON ELASTIC FOUNDATION 


FORCE = F 



Figure 1. Spherical indenter with plate on elastic foundation 



PLATE SHOWN IN TWO SECTIONS. INSIDE PORTION O-A, 
CONFORMS TO THE SPHERE WITH RADIUS Rg. 


SECTION 1 



Figure 2. Plate shown in two sections. Inside portion, O-A, 
conforms to the sphere with radius R s . 




Figure 3. Contact radius 











SECTION OF THE CIRCULAR PLATE MODEL 



i Figure 5. A section of the circular plate model 




CXI 

CO 


SIMPLIFIED FLOW DIAGRAM OF THE ITERATION PROCESS 



'Figure 7. Simplified flow diagram of the iteration process 






ELASTIC-PLASTIC ANALYSIS USING A TRIANGULAR 


RING ELEMENT IN NASTRAN 
P. C. T. Chen 

U.S. Army Armament Research and Development Command 
Benet Weapons Laboratory, LCWSL 
Watervliet, NY 12189 


SUMMARY 


An elastic-plastic triangular ring element is implemented in NASTRAN 
computer program. The plane-strain problem of partially-plastic thick-walled 
cylinder under internal pressure is solved and compared with the earlier 
finite-difference solution. A very good agreement has been reached. In order 
to demonstrate its application to more general problems, an overloaded thread 
problem for the British Standard Buttress is examined. The maximum axial and 
principal stresses are located and their values are determined as functions of 
loadings. 

INTRODUCTION ' 


The piecewise linear analysis option of the NASTRAN program provides an 
algorithm for solving nonlinear problems in material plasticity (ref. 1). 
However, the usefulness of this option is quite limited because only a few 
elements have been implemented. These include rod, tube, bar elements for 
one-dimensional problems and plate elements for two-dimensional plane stress 
problems. In a recent paper (ref. 2), the implementation of a trapezoidal 
ring element in NASTRAN for elastic-plastic analysis was described and two 
test problems of rotational symmetry were solved. The first is an infinitely 
long tube under uniform internal pressure. The NASTRAN results are in 
excellent agreement with an exact solution based on the finite-difference 
approach (ref. 3). The second problem is a thick-walled cylinder of finite 
length loaded over part of its inner surface. The NASTRAN results are in good 
agreement with another finite-element solution (ref. 4). 

In the present paper, two elastic-plastic problems of rotational symmetry 
are solved using triangular ring elements for the finite element models. The 
implementation of a triangular ring element in NASTRAN for elastic-plastic 
analysis follows the same procedures in reference 2 for a trapezoidal ring 
element. In order to test the accuracy of the implementation, the 
plane-strain problem of a partially-plastic thick-walled cylinder under 
internal pressure is solved again and compared with the finite difference 
solution (ref. 3). The second demonstrative problem in using NASTRAN 
triangular ring elements is an overloaded thread problem for the British 
Standard Buttress. A detailed discussion for this problem in the elastic 
range of loading was reported in reference 5. In the present paper, 
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uniform pressure distribution on the primary bearing surface of the thread is 
assumed and the load is applied in increments beyond the elastic limit. The 
maximum axial and principal stresses are located and their values are 
determined as functions of loadings. 

TRIANGULAR RING ELEMENTS 


The incremental displacement field employed for the triangular ring 


element are 

Au(r,z) = + 02 r + 3 3 Z > (1) 

Aw(r,z) = 04 + 05 r + 0gZ . (2) 

The transformation from grid point coordinates to generalized coordinates is 
{3} = [r pq ] {Aq} (3) 

where 

{q} T = L A u l» Aw l> Au 2 » Aw 2> ^ u 3» Aw 3 j » (4) 

{ 0 } T * L 8 l> ^ 2 , 3 3 , 04, 05 , 06 J, (5) 


and the elements of the inverse of the transformation matrix [Tfjq ] -1 are the 
coefficients of the 0 's in equations ( 1 ) and ( 2 ), evaluated at the corners of 
the element. 

The stiffness matrix is formed in the same manner as that for the 
anisotropic elastic element. The final form referred to grid coordinates is 

[k] = [r eq ] T .[K] [r 0q ] , (6) 

where 

[k] = 2ir /r[ b] t [d] [b] dzdr . (7) 

The [ b] matrix is the same as the elastic case, but now it expresses the 
incremental strains in terms of generalized coordinates 

{Ae} = [B] {0} . ( 8 ) 

The [d] matrix which relates the incremental stresses to the incremental 
strains, i.e., 

{Ao} = [d] {Ae} , (9) 

is the same as that for a trapezoidal ring element presented in reference 2 . 
In developing NASTRAN program for strain-hardening materials, we calculate 
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[ d] and obtain its inverse [d] numerically. For ideally plastic materials, 
this procedure fails and we should calculate [ d] directly using the closed 
form. 


NASTRAN IMPLEMENTATION 


The implementation of a triangular ring element in NASTRAN for 
elastic-plastic analysis follows the same procedures for a trapezoidal ring 
element (ref. 2). Changes were required in the functional modules PLA1, PLA3, 
and PLA4, which included the writing of seven new subroutines. The changes in 
PLA1 allows this module to identify the new element as a member of the 
piecewise linear element set and properly initialize the nonlinear Element 
Summary and Element Connection Property Tables. Three element stress recovery 
subroutines were added to PLA3: PSIARG, a driver for stress data recovery; 

PSRIR1 ans PSRIR2, phase I and II stress recovery routines. Element stiffness 
calculations in PLA4 require four new subroutines: PKIARG, a driver for 

nonlinear triangular ring elements in PLA4, PKRIR1 and PKRIR2, stress recovery 
routines which generate stresses for the computation of the nonlinear material 
matrix; and PKRIRS, the stiffness matrix generation routine for nonlinear 
triangular ring elements. 


PROGRAM EVALUATION 


In order to evaluate the accuracy of the computed code, the plane-strain 
problem of an elastic-plastic thick-walled cylinder under uniform internal 
pressure is solved again and compared with the finite-difference solution 
(ref. 3). The tube of outside radius 2” and inside radius 1” has been divided 
into 10 equal intervals and each interval consists of 2 triangular ring 
elements. The material constants are E = 30 x 10^ psi, v = 0.3, a Q = 1.5 x 
105 psi and the effective stress-strain curve is represented by three line 
segments connecting the four points in the (e,o) plane, (e,a/a D ) = (0.0, 0.0), 
(0.005,1.0), (0.055,1.5), (0.1, 1.5). Twenty-three load factors are chosen: 
P/o Q = 0.4323, 0.4738, 0.5125, 0.5484, 0.5818, 0.6128, 0.6415, 0.6681, 0.6925, 
0.7150, 0.7356, 0.7545, 0.7716, 0.7871, 0.8011, 0.8135, 0.8245, 0.8341, 

0.8423, 0.8493, 0.8550, 0.863, 0.87. The numerical results based on the 
NASTRAN program have been obtained. For this problem, exact solution based on 
a new finite-difference approach (ref. 3) can be used to assess the accuracy 
of the NASTRAN code. Some of the results for the displacements and stresses 
are presented graphically in Figures 1 and 2. The radial displacements at the 
inside as well as outside surface are shown in Figure 1 as functions of 
internal pressure. Figure 2 shows the distributions of radial, tangential and 
axial stress components in a partially-plastic tube when the pressure is p = 
0.7356 o 0 . As demonstrated in Figures 1 and 2, the NASTRAN results are in 
very good agreement with those based on the finite-difference approach (ref. 
3). 


192 



THREAD PROBLEM 


As an application of using NASTRAN triangular ring elements for 
two-dimensional problems, an overloaded thread problem for the British 
Standard Buttress will be examined. A finite element model is shown in Figure 
3. A detailed discussion for this problem in the elastic range of loading was 
reported in reference 5. In the present paper, uniform pressure distribution 
(P) on the primary bearing surface of the thread is assumed and the load is 
applied in ten unequal increments beyond the elastic limit. The elastic 
constants are E = 25000 Ksi, v = 0.25, and the effective stress-strain curve 
is represented by five line segments connecting the six points in the (e-o) 
plane, (e-o in Ksi) = (0.0, 0.0), (0.0048,120), (0.0122,180), (0.0167,200), 
(0.0918,210), (0.25,210). The sides A-B and C-D are constrained in the axial 
direction and the side B-C, in the axial and radial directions. 

The elastic problem is solved first and the upper limit of the elastic 
loading (p*) is found to be 54.65 Ksi. Eleven load factors are chosen as 
follows: p = 54.65, 60.59, 66.22, 71.54, 76.54, 81.23, 85.61, 89.68, 93.43, 

96.87, 100.0 Ksi. The numerical results based on the NASTRAN code have been 
obtained and the total CPU time is 59 minutes on IBM 360 model 44. The 
plastic zone at the maximum load is indicated by the shaded area in Figure 3. 
Some of the stress results are shown graphically in Figures 4 to 6. Of 
particular interest is the region along side D-E-F. The axial stress (o z ) and 
major principal stress (o^) in the boundary elements along side D-E-F are 
shown in Figures 4 and 5, for the first and last load factor, respectively. 

The maximum fillet stress (oj) occurs near E while the maximum axial stress 
always occurs at D. Finally, the maximum fillet stress and axial stress are 
plotted as functions of contact pressure in Figure 6. 


CONCLUSION 


An elastic-plastic triangular ring element has been implemented in 
NASTRAN computer program. Its accuracy has been evaluated by solving a 
simpler problem for which exact solution is available. Its application to 
more general problem has also been demonstrated by solving an overloaded 
thread problem. 
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Figure 1, Radial displacement in a pressurized tube. 



Figure 2. Stresses in a pressurized tube. 
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10 ksi 



Figure 4. Stresses in boundary elements (side D-E-F) at p = 54,65 Ksi. 
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Figure 5. Stresses in boundary elements (side D-E-F) at p = 100 Ksi. 
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DEVELOPMENT AND ANALYSIS OF THE LEARJET 54/55 FUSELAGE 
NASTRAN MODEL USING SUBSTRUCTURE TECHNIQUES 

BY 

ROBERT R. BOROUGHS, SIVAM PARAMASIVAM, AND JOANNA WERNER 


SUMMARY 


Development and analysis of the Learjet 54/55 fuselage Nastran model 
presented a real challenge considering the size of the task and the resources 
that were available at the beginning of this project. Consequently, this 
structure was broken down into several substructures to make the modeling and 
analysis effort more manageable. Since the geometry was fairly complex in some 
areas and in order to provide flexibility for future configuration studies, a 
series of local coordinate systems were used to describe the model. This work 
was accomplished using several different computer systems, the more recent of 
which were connected by a high speed data communications lines to perform the 
tasks of model generation and Nastran analysis. 


INTRODUCTION 


The Learjet Model 50 series aircraft has become the latest and largest 
member of the Learjet product line (see Figure 1). This airplane was an almost 
completely different aircraft from earlier Learjet models. The wing has been 
extended six feet from the original Learjet wing configuration, and the tip 
tanks have been removed and replaced by winglets. A paper describing the 
Nastran finite element analysis of this wing which was also used on the 28/29 
airplanes can be found in another NASA document (Ref. 1). Discussion here 
will be primarily directed toward the Nastran analysis of the fuselage and 
vertical fin. 

The fuselage on the Model 54/55 aircraft has been increased in diameter 
as well as length over previous models. Construction in the fuselage was the 
typical skin stringer arrangement with frames located at given increments to 
provide ridity to this shell structure. The windshield was similar to other 
Learjet windshields with the large unsegmented stretched acrylic panels, 
except that these panels were bigger and two small side windows were installed 
on the aft perimeter of the windshield. Attachment of the wing to the fuselage 
was accomplished through eight fittings, four on each side of the fuselage, 
which provided fail safe capability in this area. The fuselage structure was 
cut out at the wing-fuselage juncture to permit wing continuity through the 
fuselage. A keel beam was installed beneath the wing to provide a more uniform 
internal load path in this region. This keel beam structure was also continuous 
forward and aft of the wing attachment. 
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Pressure bulkhead construction consisted of either sheet metal webs 
reinforced with stringers and support beams or honeycomb panels reinforced 
with structural beams. The aircraft was powered by two Garrett Ai Research TFE 
731-3 turbofan engines mounted on each side of the aft fuselage section. The 
engine support structure consisted of a forward and aft box beam which were 
also integral to two partial bulkheads. Configuration of the vertical tail 
included five spars covered by aluminum sheet and stiffened by formed sheet 
metal ribs. 


BACKGROUND 


Initial studies on the 50 series aircraft started in 1976. Nastran was 
used as preliminary design tool to study various structural arrangements and 
determine which configuration was more advantageous. As these investigations 
matured and more parameters became defined, a transition from a conceptual 
finite element fuselage model to a final configuration model was necessary. 

The final configuration fuselage model proposed was to have a fine enough mesh 
to accurately define the internal loads distribution throughout the structure 
for all load cases, yet not be overly complex so that the lead times and costs 
would be excessive. With these guidelines in mind a model of approximately 
20,000 degrees of freedom was recommended. Development of this size model 
still presented real challenges considering the computer resources that were 
available at the beginning of this effort. There was a good deal of concern 
about the logistics of model generation, management of the resulting data base, 
and the method of analysis. The computer system initially considered to carry 
out these functions was an in-house IBM 370-145, and the method of analysis 
proposed was static substructuring. If the Nastran analysis were to become 
too large for the IBM 370-145 to handle effectively, an alternate course of 
action was planned where the Nastran model would be run on a Cyber computer at 
a service bureau. 


SUBSTRUCTURE DEFINITION 


Based on previous runs on the IBM 370-145 computer, a maximum substructure 
size of 3500 degrees of freedom was established. Each fuselage node point was 
considered to have the capability of transferring all 6 degrees of freedom. 

Grid points on the outer surface were to be located at the intersection of the 
frame and stringer members, and grid points in the fuselage interior were to 
be defined by the intersection of primary structural beams and stiffeners. 

These constraints and the basic characteristics of the structure resulted in 
dividing the Nastran model into seven substructures (see Figure 2). Before the 
modeling ever began a node and element numbering system was established for the 
entire structure. Element and grid point numbers were assigned to each sub- 
structure in a progressive manner so that one substructure could later be 
easily merged or subdivided with adjacent substructures. Node numbers were 
assigned in a sequence that attempted to minimize the matrix bandwith. Element 
numbering was keyed to the grid point numbering so that both the grid points 
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and connection members had a similar numbering sequence. 


After studying the pros and cons of various coordinate systems, a geo- 
metric definition using a series of local cylinderical coordinate systems in 
conjunction with the basic rectangular coordinate system was selected. In 
this system the grid points at each frame location on the outer surface were 
described by ; a local cylinderical coordinate system unique to that frame. Grid 
points interior to the fuselage model were defined in a rectangular coordinate 
system. There were several advantages in using this type of approach. First 
of all, this made the modeling of some of the more complex geometry much easier, 
and also simplified the determination of offsets to be used with the BAR 


elements which were to define the 


frame members. 


Consideration 


was also made 


for the incorporation of a fuselage stretch and other possible modifications 
that would impact future modeling and analysis. Thus, the use of a local 
cylinderical coordinate system for each frame would permit a fuselage stretch 
by inserting a fuselage plug and redefining the affected coordinate cards with 
little or no changes to existing grid or element cards. 


COMPUTERS USED FOR MODEL GENERATION 


Several in-house computers were used in different phases of the model 
generation. The first phase of the modeling was the generation of the connect- 
ivity cards. Simple routines were written for CONROD and CQDMEM2 elements 
(Ref. 2) which given the bay ID numbers generated each connection element. 

These routines were originally run on an IBM 1130 but all these programs were 
eventually converted to run on an IBM 370-145 or on a PDP~ 11/70. A similar 
situation occurred with some frames which were constant in section around the 
circumference. A program was used that generated all four cards necessary for 
each BAR element, complete with connectivity, cross-section area, X and Y 
moments of inertia, and stress recovery coefficients, but with all offsets 
equal to zero. These values had to be added manually, but even after adding 
these values by hand, much time was saved using these routines. 

The second phase of model generation was the section property calculation. 
This was done with two different routines. The first routine used ran on the 
IBM 370-145 and was primarily designed to calculate the section properties for 
a channel section made of bent-up sheet metal. This routine could also be 
used to calculate the properties for an angle section. A second routine that 
was used calculated the properties for any extruded section used such as 
channels, I-sections, T-sections, and angles. Special provisions for a general 
section built up of rectangles and fillets were also incorporated in this 
program. This routine ran interactively on the PDP 11/70 in-house computer. 

In some cases where a frame section was comprised of two or more bent up 
sheet sections which were fastened together, the combined section properties 
had to be computed. Using the appropriate section property program to find 
the properties of each component, these properties were then input to a trans- 
formation program which converted all the properties with respect to the 
composite centroid and combined these values. This program was set up to run 
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on the in-house IBM 370-145. 


FORWARD SECTION 


The forward section began at frame 1 and extended to frame 16 just behind 
the cabin door (see Figure 2). This portion of the model was divided into 
two substructures which were identified as the nose substructure and crew sub- 
structure. The nose substructure extended from frame 1 to frame 6 (see Figure 

3) and the crew substructure extended from frame 6 through frame 16 (see Figure 

4) . The components of the forward fuselage consisted of circumferential frames 
oriented in a fuselage station reference system, stringers located approximately 
normal to the frames, and the skins which covered this framework. Along the 
bottom centerline of the fuselage was the keelbeam which runs almost the com- 
plete length of the fuselage to give additional longitudinal stiffness and to 
provide a more continuous load path around the wing. Other components of the 
forward section included a forward pressure bulkhead, a windshield, and a cabin 
door. 


The grid points defining the outside contour for these substructures as 
well as the other substructures were located at the frame stringer intersections. 
These grid points were defined in a local cylinderi cal coordinate system - 
coincident with each frame. The local coordinate systems were defined with 
reference to a basic rectangular system where the X axis runs longitudinally 
aft, the Y axis left hand outboard, and the Z axis down. The local systems 
were oriented with R radially outward from the centerline, e counterclockwise 
from the basic Y coordinate, and positive Z oriented forward. The grid points 
were numbered with even values on the left hand side and odd values on the 
right hand side. These numbers were keyed off the local system identification 
numbers to give unique ID's for each grid point. Grid points in the fuselage 
interior utilized a similar technique. 

Stringer members in the forward section were modeled using C0NR0DS, since 
these elements had small cross-sectional areas that acted primarily in axial 
load transfer. In a few cases a torsional stiffness value was included to allow 
some torsional loading. The C0NR0D IDs directly keyed off the grid point ID's. 
The frames were modeled as BARS (Ref. 2) to carry axial and bending loads. The 
BAR orientation was defined by a v vector at the element origin pointing radially 
outward. QDMEM2 elements were used to model the skin, since these panels have 
small thicknesses and generally do not carry significant local bending loads. 

The keelbeam was modeled using C0NR0DS along the four corners with SHEAR 
elements (Ref. 2) on vertical sides. The bottom was coincident with the outer 
skin and consequently was already defined by QDMEM2 panels. 

Stretched acrylic material was used for windshield and cabin windows. 

These members were generally quite thick and were modeled using QUAD2 elements 
(Ref. 2) which had bending capability. The forward end of the pressure vessel 
was located at frame six in the form of a pressure bulkhead. This bulkhead was 
comprised of a thin web which was supported by vertical stiffeners on the forward 
side and horizontal stiffeners on the aft side. The thin web was modeled with 
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QDMEM2 elements, and the stiffeners were modeled with BAR elements. 

The cabin door structure consisted of a framework of intersecting frames 
and stringers with skin and doublers on both the inner and outer sides of 
this framework. A double set of grid points was used to define the perimeter 
of the door. The first set of grid points defined the door frame or cutout 
in the fuselage while the second set of grid points defined the edges of the 
door. Once again the frames were modeled with BARS, the stringers were 
modeled with C0NR00S or BARS depending on the depth of the member, and the 
skin panels and doublers were modeled using QDMEM2 elements. The entire 
door was connected to the fuselage by the use of RIGID (Ref. 2) elements which 
transferred the appropriate degrees of freedom. 


MID SECTION 


The fuselage mid section began at frame 16 and continued through frame 
31 just aft of the wing trailing edge (see Figure 2). This section was divided 
into two substructures at frame 24 which was just in front of the wing leading 
edge. The substructure between frames 16 and 24 was identified as the cabin 
substructure (see Figure 5), and the substructure between frames 24 and 31 
was called the center substructure (see Figure 6). The shell structure 
for these sections consisted of frame-stringer-skin type construction and 
contained the keel beam, the escape/baggage door, passenger windows, the 
frame 24 partial bulkhead, a pressurized baggage floor, the aft pressure 
bulkhead, and a portion of the fuselage fuel cell bay. 

The grid point locations and modeling of the fuselage shell structure 
for the mid section was accomplished in basically the same manner as described 
in forward section modeling discussion. The keel beam which runs almost the 
entire length of the fuselage was an open box type structure forward of frame 
24, but aft of frame 24 and the wing leading edge region the keel beam has 
become a completely closed box. Modeling of the keel beam in this region was 
accomplished using two dimensional Nastran elements. CONROD members were used 
to model both the upper caps, lower caps, and vertical stiffeners. SHEARS 
were used to model the vertical webs while QDMEM2 elements were used to model 
the horizontal webs. 

Window cutouts were reinforced with doublers as well as frame members in 
certain locations. The doubler panels were modeled using QDMEM2 elements and 
were connected to the same grid points as the outer skin panels. The windows 
were made of the same stretched acrylic material as the windshield but only 
thinner, and these panels were also modeled with QDMEM2 elements. 

Between frames 22 and 25 and stringers 6 through 15 on the right side of 
the fuselage was located the escape/baggage door. A double set of grid points 
was used to define the perimeter of the escape/baggage door similar to the 
modeling performed on the main cabin door. Apart from the two close out 
frames at each end of the door, there were also two inner frames which were 
adjacent to the window in the door. The door cutout in the fuselage was 
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reinforced with doublers and additional frame members and intercostal s. The 
door frame and intercostal members were modeled using BAR elements, and the 
door skin was modeled using QDMEM2 elements. Connection of the escape/baggage 
door to the fuselage was achieved by using RIGID (Ref, 2) elements between the 
appropriate degrees of freedom at the door hinge and pinned attachments. 

The partial bulkhead at frame 24 served as a close out for the wing cutout 
in the fuselage. This structure consisted of a thin web supported by hori- j 
zontal and vertical stiffeners. QDMEM2 elements were used to model the web, 
and BAR elements were used to model the stiffeners. The top edge of this 
bulkhead provided the forward support point for the pressurized baggage floor 
which was made of sandwiched honeycomb plate. Due to the thickness and 
bending characteristics of the honeycomb, this structure was represented in 
the Nastran model by QUAD1 elements (Ref. 2). The baggage floor was attached 
on both sides to a longeron which ran the length of the wing cutout in the 
fuselage. These longerons also helped to close out the wing cutout and pro- 
vided a redistribution path for the internal loads. Since the longeron was a 
large member, BAR elements were used to similate this structure in the model 
instead of CONRODS. 

Attachment of the wing to the fuselage was accomplished through four 
fittings on each side of the fuselage in the wing cutout region. The frames 
located through this section were all double frames. This was done to provide 
increased stiffness and an adequate load path for the wing reactions. These 
double frames were all modeled with a single ring of BAR elements to conserve 
degrees of freedom. Since the wing structure was not simulated in this model, 
the model was constrained at each of these wing attachment points. The frame 
28 structure also provided the support point for the aft pressure bulkhead. 

This bulkhead also served as the forward retainer for the fuselage fuel cell, 
and was constructed of sandwiched honeycomb plate. QUAD1 elements were used 
to model this bulkhead as was previously done for the pressurized baggage floor. 


AFT SECTION 


The aft fuselage section began just behind the trailing edge of the wing 
and was divided into two substructures (see Figure 2). These substructures 
were referred to as the fuel cell (see Figure 7) and tailcone substructures 
(see Figure 8) and had a mutual boundary at frame 39 which was between the 
two aft baggage doors. Some of the basic assemblies of the fuel cell section 
included the fuel cell support structure and aft fuel cell bulkhead, the 
engine support structure, and the aft keel beam structure. Major features 
in the tailcone substructure consisted of the aft baggage compartment and the J 
vertical tail support structure. Frames in both substructures were generally 
oriented in a vertical position with the exception of the frames that attached 
to the vertical fin which were oriented parallel to the spars in the vertical 
fin. 


Nastran elements used to model the skin, stringer, and frame elements 
were the same as those used in the other substructures. There was one full 
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bulkhead in the aft section and several partial bulkheads. The full bulkhead 
served as the rear boundary of the fuselage fuel bay and was located just in the 
front of the forward engine beam,. This bulkhead was constructed of flat plate 
reinforced by structural beams. QUAD2 elements (Ref. 2} were used to model 
the plate structure, and BAR elements were used to model the structural rein- 
forcement beams. One partial bulkhead was located at the aft end of the wing 
fuselage cutout at frame 31 and served as a close out as well as a fuel 
retainer for the fuselage fuel bay. Construction of this partial bulkhead 
was also a flat panel with structural support beams, but the reinforcement 
beams were generally smaller and there were more of these members than in the 
aft fuel bulkhead. The webs for this partial bulkhead were modeled with QUAD2 
elements, and the support beams were modeled with BARS. 

Since grid points in the bulkhead mesh did not always match the grid 
points on the fuselage outside contour in an even manner, the interface of 
cross support beams on the bulkheads with the fuselage frames was not always 
easy, to simulate in the Nastran model. MPC equations were used initially to 
relate the displacement of the grid points in the bulkhead which were adjacent 
to the outside fuselage contour to the displacement of the grid points on the 
outside contour. The results of the first Nastran debugging runs with these 
MPC equations revealed the reactions did not satisfy equilibrium conditions. 
However, when the MPC equations were removed from the model , then equilibrium 
was satisfied. This problem was reported to COSMIC, and a sample Nastran run 
with the MPC equations and faulty equilibrium conditions was sent to the 
COSMIC staff. After reviewing the Nastran run, COSMIC indicated that the 
problem appeared to be related to writing an MPC equation for grid points 
which were in different coordinate systems whose degrees of freedom were not 
on parallel reference axes. The coordinate system used for the grid points 
on the fuselage outside contour was a local cylinderical coordinate system 
with the origin at the maximum breadth line of the fuselage while the 
coordinate system used for the grid points on the fuselage interior was a 
rectangular coordinate system. Consequently, MPC equations were discarded for 
this application and tailored BAR elements were used instead. 

Two of the other partial bulkheads were located at the forward and aft 
engine beam supports. These engine beams were box beams that were continuous 
through the fuselage from right to left. Both beams were curved so that the 
outboard ends, or engine attach points, were higher than that portion of the 
engine beam on the aircraft center Tine. A double frame was also installed 
at each of these locations to provide. greater support for the engine loads. 

In order to simplify the modeling and conserve degrees of freedom, the 
combined section properties of the double frame were calculated and one ring 
of BAR elements was used to model both frames. The partial bulkheads in 
this region extended from the engine beam on the bottom to the double frame on 
the sides and on the top. BAR elements were used to model the engine beams, 
and QUAD2 members were used to model the web panels in the partial bulkhead. 
There were only a few reinforcement members in the partial bulkhead, and these 
elements were modeled using BAR elements. 

The intersection of the engine beams and the stringer members on the 
fuselage contour were such that the grid points for these elements did not 
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coincide at the true intersection point. Since there was internal load transfer 
capability of all six degrees of freedom at this point, a method of relating 
this interaction was necessary, RIGID elements (Ref, 2) were used to connect 
these degrees of freedom, and the results appeared to be satisfactory. 


VERTICAL FIN 


The seventh substructure in the fuselage model was the vertical fin installa- 
tion which also included a model of the rudder (See Figure 9). A model of 
the horizontal stabilizer was constructed, but was not included in this data 
base in order to minimize the degrees of freedom and cut down on computer run 
costs. Horizontal stabilizer loads were applied to the vertical fin at the 
attach points where the horizontal stabilizer and vertical fin joined. The 
vertical fin consisted of five spars and eight ribs covered with aluminum 
sheet while the rudder contained one spar and eleven ribs with an aluminum 

covering. The vertical fin was joined to the tailcone at each of the five 

spars as well as through attachments on the skin. Spar connections were made 
at frames 43, 45, 46, 47 and 48 which were canted to accommodate this interface. 

A separate local rectangular coordinate system was adopted for the grid 
locations for both the vertical fin and rudder to facilitate the modeling. 

The spar caps in the region of the fuselage frames were modeled using BAR 
elements while the remainder of the spar caps and rib caps were modeled using 
CONROD elements. Skin panels were modeled using the QDMEM2 member, and SHEAR 
elements were used to represent the spar and rib webs. Attachment of the 
rudder to the vertical fin was accomplished through two hinges which were 
modeled using CONRODS. The torque load in the rudder was restrained by a 

torque tube attached to the bottom of the rudder which in turn attached to 

fuselage frame 48 and was modeled using BAR elements. The spar caps and rib 
caps were again modeled using CONROD elements while the skin was modeled using 
QDMEM2 panels. CSHEAR elements were used to model the rudder leading edge 
ribs while TRMEM elements were used for the ribs between the hinge line and 
trailing edge of the rudder. 


GRAPHICS SYSTEMS 


Plots for each of the completed substructures were obtained using an 
in-house plot package developed for Nastran on an in-house PDP 11/70 mini 
computer. Since this program was originally written for plots with grid 
points in the basic coordinate system, an expansion to the software was made 
to incorporate multiple local coordinate systems. The user could access this 
routine through a Tektronix CRT, and the substructure would be displayed in 
a three-dimensional view on the screen. Different views could be selected by 
specifying the appropriate rotational angles about each of the three principal 
axes. Hard copy plots of very fine resolution then could be obtained by 
spooling the desired view to a Versatec plotter. This approach allowed the 
user to correct any noticeable geometry or connectivity errors before going 
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to the initial Nastran finite element analysis. 


LOADS 


Several different load conditions and combinations of load conditions were 
applied to the model. The basic types of loads used during this analysis were 
pressure loads, down bending loads, up bending loads, and side bending loads. 

An ultimate and limit internal pressure was applied separately as well as in 
conjunction with the ultimate down bending and ultimate up bending load cases. 
Side bending loads were applied to the vertical fin and the rudder. 


COMPUTERS USED FOR ANALYSIS 


• 

Further error correction was achieved by running the substructures with 
the Nastran software on the IBM 370-145. However, after the Nastran debugging 
process was completed, the first few substructures ran so long that even over 
night turn around became a problem. This difficulty appeared to be signifi- 
cantly influenced by increased usage of the computer not only by structures 
personnel but also by other departments. One of the alternatives considered 
as a solution for this situation was to break the model into smaller sub- 
structures and increase the number of substructures from seven to eleven sub- 
structures. A second alternative was to run the model on a more powerful out- 
side computer. At that point in time the second of the two alternatives 
seemed to be more satisfactory in order to meet schedule requirements. Several 
different outside sources were examined, and the selection process finally 
narrowed down to Control Data's Cybernet System running Nastran on a Cyber 175 
computer. 

The first complete computer runs for the 54/55 Nastran fuselage analysis 
were made successfully on the Cybernet System using Nastran Superelement 
analysis. In order to cut down on the transmission of data over the phone' 
lines, a data base was set up for the fuselage model at the computer site. 

This data base permitted the running of several different load cases at 
reduced cost by transmitting only the JCL and loads data between the Lear jet 
terminal and the Cyber 175. 


EXPERIMENTAL RESULTS 


Only a portion of the static test ’..program had. been completed, and the 
correlation pf experimental data;. with analytical" resul ts had been under way 
a short time when thisTpaper had to be submitted to the. publ isher: i. /Consequently , 
just a- few figures were available showing the comparison between' the Nastran 
results and the static test data. A sample of two load conditions have been 
shown on the following .pages. These two load conditions were the vertical fin 
side bending load case and the limit pressure loading. Plots for the side 
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bending load case have been shown in figures 10 and 11, and the correlation 
between the Nastran values and experimental data for the limit pressure load 
case have been shown in figure 12. 


CONCLUDING REMARKS 


Modeling of the Learjet 54/55 fuselage was based upon use of a multiple ! 
local coordinate system which proved to have significant advantages and 
worked very satisfactorily. The fuselage was analyzed using Superelement analysis 
where the structure was divided into seven superelements. This approach made 
the analysis of this size structure much more manageable and easier to perform. 
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FIGURE 3 - NOSE SUBSTRUCTURE 



214 






FIGURE 6 - CENTER SUBSTRUCTURE 
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FIGURE 7 - FUELCELL SUBSTRUCTURE 
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FIGURE 8 - TAILCONE SUBSTRUCTURE 
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FIGURE 9 - VERTICAL FIN SUBSTRUCTURE 
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COMPARISON OF FINITE ELEMENT ANALYSES OF A 


PIPING TEE USING NASTRAN AND CORTES/SA 

Antonio J. Quezon and Gordon C. Everstine 
David W. Taylor Naval Ship Research and Development Center 

SUMMARY 


A comparison of finite element analyses of a 24" x 24" x 10" piping tee 
was made using NASTRAN and CORTES/SA, a modified version of SAP3 having a 
special purpose input processor for generating geometries for a wide variety of 
tee joints. Four finite element models were subjected to force, moment, and 
pressure loadings. Flexibility factors and principal stresses were computed 
for each model and compared with results obtained experimentally by Combustion 
Engineering, Inc. Of the four models generated, the first was generated from 
actual measured geometry using GPRIME, a geometric and finite element modeling 
system developed at DTNSRDC. The other three models were generated from an 
idealized tee using the data generator contained in CORTES/SA. 

The generation of an idealized tee proved to be very easy and inexpensive 
compared to generation from actual geometry, and, when analyzed by NASTRAN, 
proved adequate. Results from the NASTRAN analyses were in good agreement 
with experimental results for all loadings except internal pressure. The 
CORTES/SA analyses gave good results for the internal pressure loading, but 
poorer results for out-of-plane bending moments or forces resulting in out-of- 
plane bending. Two of the basic load cases in CORTES/SA were found to contain 
errors that could not be easily corrected. A cost comparison of NASTRAN and 
CORTES/SA showed NASTRAN to be less expensive to run than CORTES/SA for 
identical meshes. Overall, considering modeling effort, cost, and accuracy, it 
is concluded that tees can be easily and accurately analyzed by NASTRAN using 
an idealized mesh generated by CORTES/SA. 


BACKGROUND 


The designer of a piping system requires a knowledge of the deflections 
and stresses caused throughout the system by anticipated service loads. Of 
particular interest are critical components such as elbows and tees. 

The purpose of this paper is to summarize the results of a study (ref . 1) 
made recently to assess the effectiveness of the finite element method (FEM) in 
predicting flexibility factors and stresses in piping tees subjected to force, 
moment, and pressure loadings. A similar study (ref. 2), performed for piping 
elbows, indicated that very good agreement could be expected between FEM 
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analysis and experiment. Tees, although conceptually no more difficult to 
analyze than elbows, are considerably more complicated geometrically. A 
reducing tee, for example, has in the crotch region a fillet with a variable 
radius of curvature as well as variable thickness. Moreover, the adjacent 
straight sections may not be cylindrical. Thus, geometrical idealizations of 
tees, although plausible, may be incorrect. 

The finite element analyses described here involve idealized models as well 
as a model based on actual measured geometry. Two computer programs were used 
for the analyses: NASTRAN, a widely-used general purpose finite element 

structural analysis program, and CORTES/SA (ref. 3), a special purpose finite 
element tee analysis program based on SAP3 and written at the University of 
California at Berkeley under the sponsorship of the Oak Ridge National 
Laboratory.* This series of analyses was designed to provide information on 
the sensitivity of the results to various mesh densities as well as on the 
adequacy of the assumed idealizations. 

In this paper, the program CORTES/SA will be referred to by the 
abbreviated name "CORTES". 


STATEMENT OF THE PROBLEM 


Combustion Engineering, Inc., performed an experimental stress analysis 
(ref. 4) on an ANSI B16.9 carbon steel tee designated T-12. Pipe extensions 
were welded to the branch and run ends of the tee, and the resulting assembly 
was placed in a load frame. One of the run ends was built in to represent a 
fixed end, and the other run end and the branch end were used to apply six 
orthogonal moments and five orthogonal forces. Internal pressure was also 
applied. Table 1 and Figure 1 summarize the applied loads. Load case 4 (F3X) 
was not tested because of strength limitations of the load frame. Stress data 
for all twelve load cases were gathered from strain gages fixed on specific 
rows on the tee (Figure 2) and were plotted against normalized surface distance. 

The tee analyzed was a reducing tee with a 24- inch-diameter run end and a 
10.75-inch-diameter branch. Loads to the run were applied at the free end of 
the run pipe extension, 173 inches from the branch-run intersection (Figure 1) . 
Loads to the branch were applied at the end of the branch pipe extension, 

77 inches from the branch- run intersection. The run pipe extension consisted 
of 24-inch-diameter schedule 40 (0.687-inch nominal wall thickness) carbon 
steel piping. The branch pipe extension consisted of 10 . 75- inch-diameter 
schedule 40 (0.365- inch nominal wall thickness) carbon steel piping. 

The finite element analyses of the tee simulated these loading conditions 
so that stresses at selected locations could be compared to the experimental 


The CORTES package of computer programs is distributed as program number 759 
by the National Energy Software Center (NESC) , Argonne National Laboratory, 
9700 S. Cass Avenue, Argonne, Illinois 60439. 
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results. For most load cases, the strain gage rows (Figure 2) selected for 
comparison were those on which the peak stresses occurred. 


ANALYSES PERFORMED 


NASTRAN analyses were performed for the first three models generated, and 
CORTES analyses were performed for the third and fourth models. These five 
finite element analyses are summarized in Table 2. In the abbreviations Nl, 

N2, N3, C3, and C4 used to identify the analyses, the first character (N or C) 
indicates the analysis program used (NASTRAN or CORTES), and the second 
character indicates the mesh used. A typical mesh generated by CORTES is 
shown in Figure 3. In general, a higher mesh number corresponds to a finer 
mesh, either overall or in selected key regions of the tee. 

. The NASTRAN analysis of Mesh 1 was the only analysis performed for a model 
generated from actual measured geometry. The remaining analyses were performed 
either by NASTRAN or by CORTES on meshes generated by CORTES assuming an 
idealized geometry. In all cases only one-fourth of the actual tee was modeled 
due to symmetry. 

For the NASTRAN analyses, the tee, including pipe extensions, was modeled 
with plate (NASTRAN QUAD2) elements. Flexible beam (BAR) elements were 
arranged in a spoke formation radiating from an imaginary point in the center 
of the cross section at the ends of the tee branch and run to facilitate the 
calculation of the average rotation of these cross sections. Rigid (RIGD1) 
elements were defined at the ends of the pipe extensions for use in load 
application. The loads were applied to a point in the center of the rigid 
cross section at the ends of the pipe extensions. 

In the CORTES analyses, the tee and pipe extensions were modeled using an 
8-node hexahedral element. This element, designated ZIB8R9, is a modification 
of the standard Zienkiewicz-Irons isoparametric element and, according to 
Gantayat and Powell (ref. 3), has bending properties superior to those of the 
unmodified isoparametric element . 

Mesh 1 was modeled from actual geometry as specified in the Combustion 
Engineering, Inc., report (ref. 4) which tabulated coordinates of points on the 
outer surface of the tee and thicknesses at these points. From these 
digitized data, a general B-spline surface was fitted through the supplied 
points using the geometric and finite element modeling processor GPRIME 
(refs. 5 and 6). Once this geometric model was defined, GPRIME was used to 
generate a finite element mesh which included the effects of variable 
thicknesses. 

Meshes 2 through 4 were modeled as idealized tee joints using the auto- 
matic mesh generation routine in CORTES. The tee joint is idealized by 
shallow cones representing the branch and run portions of the tee, connected 
to each other through an analytically defined transition fillet (ref. 3). 
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STRESS RESULTS 


The results of primary interest are normalized principal stress values 
for elements in particular locations on the tee. The peak normalized 
principal stress was plotted against surface distance ratio for each load case 
and compared to the experimental results obtained by Combustion Engineering, 

Inc. (ref. 4). 

The tee analyzed by Combustion Engineering, Inc., was heavily instru- 
mented, both internally and externally, with strain gages in two of the four 
quadrants. The gages in each quadrant were arranged in six rows as shown in 
Figure 2. Since the peak stresses for most load cases occurred on row 1 or 
row 6, analytical and experimental results were compared for these rows only. 

For each load case, the analytical results for principal stresses were 
normalized by a stress calculated from beam theory, as indicated in Table 1. 

The normalized principal stresses were then plotted against the surface distance 
ratios of the elements lying on row 1 and row 6. 

Stress plots for several typical load cases are shown in Figures 4 
through 8. (Ref. 1 contains plots for all load cases.) All finite element 
curves are smoothed slightly by fitting B-spline curves (refs. 7 and 8) through 
the computed values, which are located at element centroids for the NASTRAN 
results and at grid points for the CORTES results. 

FLEXIBILITY FACTORS 


Two ambiguities were encountered in comparing computed flexibility factors 
with experimental results obtained by Combustion Engineering, Inc. These 
ambiguities involved the definition of flexibility factors and the way in which 
the rotation of branch or run end cross sections was measured. Combustion 
Engineering, Inc., defined the flexibility factors as 


k = 


0 - 0 
meas corr 


6 

nom 


( 1 ) 


where 


0 

meas 


measured rotation at an intermediate location on the pipe 
extension 


0 

corr 


8 

nom 


rotation correction computed by simple beam theory for the length 
of pipe between the tee weld line and the location at which the 
rotation is actually measured 

nominal rotation computed by simple beam theory for the distance 
between the tee weld lines where 
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(for bending moments) (2) 

(for torsional moments) (3) 

(for point loads) (4) 

Since Combustion Engineering, Inc., could not measure the actual rotations 
at the branch and run end cross sections, measurements were made at other 
locations on the pipe extensions and then corrected to the branch and run ends 
using simple beam theory. On the other hand, the NASTRAN analyses used very 
flexible beam elements radiating from an imaginary point in the center of the 
branch and run end cross sections to the points on the circumference of the 
branch and run ends. This modeling technique allowed an approximate average 
rotation for the cross sections to be easily obtained for the imaginary center 
point. However, because plane sections do not, under loading, remain plane, 
there is no single rotation for a section, so that different methods for 
computing rotations will yield different results. 

For the computation of flexibility factors from the NASTRAN results, the 
relation 

k = ~~ ' (5) 

nom 

was used, where 


e 

nom El 

e .= S*. 

nom JG 
PL 2 

8 = 
nom 2EI 



0 

nom 


= computed relative rotation of end "a" with respect to end "b" 

= nominal rotation computed by beam theory for the rotation of end 
"a" with respect to end "b" 


Flexibility factors were computed for the free branch and run ends with 
respect to the fixed run end for each load case except for F2X (an axial load 
on the run) and internal pressure, neither of which causes any significant 
rotation. For example, the flexibility factor for a rotation about the X-axis 
of the branch end with respect to the fixed run end is denoted by kx 3 i> where 
the X in the subscript represents the axis of rotation, the 3 represents the 
branch end, and the 1 represents the fixed run end. For each load case, 
flexibility factors for each cross section were computed. 


Table 3 compares the flexibility factors computed from the three NASTRAN 
analyses to the experimental values. The computed flexibility factors compare 
reasonably well for most load cases, an exception being k^2i for load case- 5 
(F3Y) of N2. Combustion Engineering, Inc., did not compute flexibility factors 
for this load case because the stresses and deflections were considered too 
small to give reliable answers. The displacements computed in the three 
NASTRAN analyses for load case 5, however, did not appear to be significantly 
smaller than those of the other load cases, although the run end of the tee did 
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warp severely in all three analyses. Since the distortions in all three 
analyses were similar, it appears to have been due to chance that the flexibil- 
ity factors for N1 and N3 were not also negative for this load case. This 
implies that any method used to compute a single rotation of the run end is 
inadequate for severely distorted cross sections. Moreover, the usefulness of 
a flexibility factor when severe cross-sectional distortion occurs is 
questionable. 

In general, a negative flexibility factor, whether arising from experiment 
or analysis, is physically impossible, since such a factor implies a rotation 
in a direction opposite to that of the applied moment. Negative values can 
arise experimentally whenever rotations measured at one location have to be 
"corrected" (using beam theory) to yield rotations elsewhere. Negative values 
can result from a finite element analysis whenever severe cross-sectional 
distortion occurs, in which case the usefulness of an "average" rotation of the 
cross section is in doubt. 


DISCUSSION OF RESULTS AND CONCLUSIONS 


The three NASTRAN analyses of the tee joint were generally in very good 
agreement with the experimental results and accurately predicted peak stresses 
for most loadings except load cases 3 (M3Z) and 13 (pressure). Also, as 
expected, the agreement with the experimental results improved with finer 
meshes. In general, the two CORTES analyses were slightly less accurate than 
the NASTRAN analyses except for load case 13 (pressure) . We are unable to 
explain this behavior. While the CORTES results for pressure loading were 
significantly better than NASTRAN's, the results for load cases 1 (M3X) , 

8 (M2Y) , 10 (F2X) , and 12 (F2Z) were worse. Note that most of these load cases 
involve either out-of-plane bending moments or forces resulting in out-of-plane 
bending. The CORTES analyses of load cases 6 (F3Z) and 7 (M2X) were also found 
to contain errors in formulation and coding which could not be easily 
corrected. 

The preparation of the NASTRAN model of mesh 1 (called Nl) was the most 
time-consuming and expensive of all the models, since this mesh was generated 
from actual geometry. Although the Nl calculations for all load cases except 
pressure (Figure 8) are in very good agreement with the experimental results, 
they are not significantly better than those obtained from the other analyses, 
so the extra effort is not justified. 

In a comparison of NASTRAN and CORTES analyses of an identical mesh 
(Mesh 3), the NASTRAN results (N3) were more consistent and predicted peak 
stresses more accurately than CORTES for ten of the twelve load cases . Only 
for M3Z and pressure (Fig. 8) did C3 do better than N3. Although N3 was less 
expensive than C3 in computer costs, it required slightly more time for input 
preparation. 

Since N3 was in generally better agreement with experimental results than 
C3, a coarser mesh (Mesh 2) was also analyzed by NASTRAN (N2) and compared 
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with C3. In all but three of the load cases, M3Z, F3Y, and pressure, N2 was 
again in better agreement with experimental results than C3. Computer costs 
from N2 were significantly less than those for C3, as indicated in Table 2. 

In an effort to obtain better results from CORTES, a much finer mesh 
(Mesh 4) was generated and analyzed, so that the results could be compared with 
N3. This time, overall performance was about equal for the two analyses, 
although C4 achieved better results than N3 for M3Z, F3Y, M2Z, F2Y, and 
pressure. 

In conclusion, it is apparent that GPRIME, although well-suited in general 
to the generation of tee meshes based on actual geometry, is more difficult 
and time-consuming to use than the special purpose idealized tee generator 
contained in CORTES. Models based on actual geometry also require geometric 
data that would probably not be generally available to the analyst. For these 
reasons, CORTES generation of a finite element model based on idealized 
geometry appears to be acceptable. However, if an analyst is interested in an 
F3Z or an M2X loading, CORTES should not be used as the analyzer because the 
program currently contains errors in the coding of these two load cases. Also, 
as shown by the comparison of N3 with C4, CORTES requires a mesh about 20% 
finer to obtain results as accurate as NASTRAN. 

Overall, considering modeling effort, cost, and accuracy, it is concluded 
that tees can be easily and accurately analyzed by NASTRAN using an idealized 
mesh generated by CORTES /SA. 
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TABLE 1 - SUMMARY OF APPLIED AND NORMALIZED LOADS 


Load 

Case 

Applied 

Load 

Nominal 

Stress 

Normalized 

Load 

1. M3X 

4.29E5 in- lb 

M3X 

Z b 

29.91 

2. M3Y 

-6.03E5 in- lb 

M3Y 

•7 

% 

29.91 

3. M3Z 

5.98E5 in-lb 

M3Z 

Z b 

29.91 

5. F3Y 

4.0E4 lb 

F3Y 

11.91 

6. F3Z 

5.58E3 lb 

77F3Z 

3.884E-1 

7. M2X 

4.9E6 in- lb 

M2X 

Z 

r 

285.0 

8. M2Y 

-7.54E6 in- lb 

M2Y 

Z 

r 

285.0 

9. M2Z 

3.40E6 in-lb 

M2Z 

Z 

r 

285.0 

10. F2X 

-6.28E5 lb 

F2X 

A 

r 

50.3 

11. F2Y 

2.01E4 lb 

173F2Y 

Z 

r 

1.6474 

12. F2Z 

2.46E4 lb 

173F2Z 

Z 

r 

1.6474 

13. P 

600 psi 

/ PD ° \ 

5.725E-2 

\ 2t ) x 


Notes; 

1. Load case 4 (F3X) was not tested. 

2. The "normalized load" is computed by dividing the experimentally 
applied load (column 2) by the nominal stress (column 3) . 

3. Subscripts "r" and "b" above denote "run" and "branch", respectively. 
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TABLE 2 - COMPARISON OF FINITE ELEMENT ANALYSES: 
NAS TRAN vs. CORTES 



N1 

N2 

N3 

C3 

C4 

NASTRAN or 

CORTES 

Analysis 

NASTRAN 

NASTRAN 

NASTRAN 

CORTES 

CORTES 

Idealized or 

Actual 

Geometry 

Actual 

Idealized 

Idealized 

Idealized 

Idealized 

Number of 
Elements 

432 

420 

525 

549 

626 

Number of 
Nodes 

484 

473 

583 

609 

689 

Number of 
Degrees of 
Freedom 

2525 

2462 

3047 

3458 

3958 

Total CP 
Seconds 
(CDC 6400) 

2213 

2200 

3135 

3310 

4748 

Cost 

$228 

$226 

$335 

$421 

$605 
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TABLE 3 - COMPARISON OF FLEXIBILITY FACTORS OF 
EXPERIMENTAL RESULTS TO NASTRAN RESULTS 


Load 

Case 

k 

Subscript 

Experiment 

N1 

N2 

N3 

1. M3X 

X21 

-0.8 

0.82 

1.00 

0.85 

X31 

1.8 

2.40 

4.00 

2.73 

2. M3Y 

Y21 

0.5 

0.72 


0.77 

Y31 

-0.3 

0.32 

0.33 

0.32 

3. M3Z 

Z21 

0.5 

0.73 

1.03 

0.93 

Z31 

0.9 

0.90 

0.85 

0.84 

5. F3Y 

Z21 


1.22 

-1.35 

1.53 

Z31 


1.97 

0.51 

2,08 

i 

6. F3Z 

X21 


0.85 

1.09 

0.88 

X31 

1.8 

2.97 

4.94 

3.40 

7. M2X 

X21 

-0.4 

0.82 

1.00 

0.85 

X31 

-0.5 

0.82 

1.00 

0.85 

- 

8. M2Y 

Y21 

0.7 

0.72 

0.76 

0.77 

Y31 

0.6 

0.72 

0.76 

0.77 

9. M2Z 

Z21 

0.9 

0.73 

1.01 

0.93 

Z31 

0.8 

0.73 

1.01 

0.93 

11. F2Y 

Z21 

0.8 

0.73 

1.01 

0.93 

Z31 

0.8 

0.83 

1.08 

1.01 

12. F2Z 

Y21 

0.7 

0.72 

0.76 

0.77 

Y31 

0.7 

0.91 

0.89 

0.94 
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Source: Figure 6, Ref, 4 


Figure 1 - Schematic View of Applied Loads for Test Tee 





Source: Figure 4, Ref. 4 


Figure 2 - Location of Strain Gage Rows on Test Tee 
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Figure 3 - Finite Element Model of Tee Using Idealized Geometry; Mesh 2 



Surface Distance Ratio 


Figure 4 - Normalized Stress Intensity for Load Case 1 (M3X) , Row 1 
Major Principal Stress on Outer Surface 





Figure 5 - Normalized Stress Intensity for Load Case 7 (M2X) , Row 1, 
Major Principal Stress on Outer Surface 
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Stress 
Intensity , 


Nl: NASTRAN, Mesh 1 
N2: NASTRAN, Mesh 2 
N3: NASTRAN, Mesh 3 
C3 : CORTES, Mesh 3 
C4: CORTES, Mesh 4 


m 





-31 1 

-1 0 1 

Surface Distance Ratio 


Figure 6 - Normalized Stress Intensity for Load Case 9 (M2Z) , Row 1, 
Minor Principal Stress. on Inner Surface 
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Nls NASTRAN* Mesh 1 

N2: NASTRAN, Mesh 2 

N3: NASTRAN, Mesh 3 

C3: CORTES, Mesh 3 

C4: CORTES, Mesh 4 
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